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Kxperimental Investigation of Turbulent 
Boundary Lavers in “Hypersonic Flow 


R. KENNETH LOBB,?7 


EVA M. WINKLER,f ann JEROME 


PERSHt 


1.8. Naval Ordnance Laboratory, Silver Spring, Maryland 


ABSTRACT 


investigation in the NOL 12- by 12-cm 
4 on naturally turbulent boundary 
and 7.7 with and without 


The Reynolds 


The results of an 
Hypersonic Wind Tunnel No 
lavers at Mach Numbers of 5.0, 6.8, 
steady-state heat transfer to the wall a 
Number 
varied from 


re given 


boundary-laver momentum thickness was 


Measurements of pitot 
and rates of heat transfer 


based on 
5,000 to 13,000 and static 
pressures, total and wall temperatures, 
made it possible to compute velocity profiles, temperature pro 
and boundary-layer parameters without resorting to any 
boundary-layer ve 


files, 
The turbulent portion of the 
was found to differ from the 


assumptions 


locity profile incompressible flow 
velocity profile in the laminar sublayer 


Local skin friction coefficients as calcu 


logarithmic law. The 
has a linear distribution 
lated from the velocity gradients at the wall and also from the 
heat-transfer coefficients using a modified Reynolds analogy are 
in accord with the results of other experimenters at lower Mach 


Numbers for the same Reynolds Number 


INTRODUCTION 


7. MAJOR PART of the drag of slender missiles at 
hypersonic** Mach Numbers is contributed by 
skin friction rather than the pressure distribution on the 
body. Furthermore, one of the limiting factors in 
hypersonic missile design may well be the extreme sur- 
face temperatures produced by the frictional heating of 
the missile. Knowledge of the boundary-layer proper- 
ties at high Mach Numbers is therefore essential to the 
missile designer because friction drag and heat transfer 
can then be determined. 
Aerodynamics Session, Twenty-Second An 
nual Meeting, IAS, New York, January 25-29, 1954 

* This investigation was sponsored by the U.S. Navy, 
The assistance of R. Garren during the tests is grate 


Presented at the 


Bureau of 


Ordnance 
fully acknowledged 

+ Chief, Hypersonics Branch 

: Physicist. 
defined as flow at Mach 


** Hypersonic is herein arbitrarily 


Numbers in excess of 5.0 


Despite the major significance of these problems, rela 
tively little experimental evidence in support of the 
various theoretical treatments is available. For the 
case of laminar boundary-layer flows, the lack of experi 
mental support for the theories is not generally regarded 
as a major inadequacy. This is because laminar bound 
ary-layer theory is regarded as more or less exact. In 
sharp contrast, however, most of the analyses of com 
pressible turbulent boundary layer flow are, in some 
manner, necessarily based on experimental results ob 
tained in wholly incompressible flow. In view of this, it 
is not surprising that the uncertainties in such analyses 
have led to large discrepancies in the prediction of skin 
friction and heat transfer. Chapman and Kester’ have 
pointed out that even if we restrict ourselves to the 
several turbulent skin friction theories that agree with 
each other and with the existing experimental data to 
within 5 per cent, they still differ greatly at hypersonic 
Mach Numbers. 

The 


compressible 


contains several papers concerning 
boundary-layer skin 


among these are the 


literature 
turbulent 
Noteworthy 


friction 
measurements. ' 
direct force measurements of Coles* and Chapman and 
Kester.' All of these investigations were primarily in- 
terested in skin friction data alone, which of course is 
one of the end products desired. However, they provide 
little information about the fluid properties within the 
turbulent boundary layer, from which a semiempirical 
theory may possibly be derived. The measurements 
have in general been made on insulated walls and in the 
Mach Number range up to 4.5 

The present paper contains detailed experimental 1n- 
formation on turbulent boundary layers in hypersonic 
The purpose of this investigation was to extend 
also to 


flow. 
the Mach Number range of available data and 
attempt to provide a deeper insight into the character 


1544808 
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istics of a turbulent boundary layer in compressible 


flows. 
SYMBOLS 
cf = local skin friction coefficient based on free-stream con- 
ditions 2ty/p oto? 
cfi_ = incompressible local skin friction coefficient for zero 
heat transfer based on free-stream conditions 
Cp = specific heat at constant pressure 
k = thermal conductivity 
MI = Mach Number 
Nu = Nusselt Number for stagnation temperature probe" 
Pr = Prandtl Number 
py) = stagnation pressure 
p)’ = pitot pressure 
p = static pressure 
Re = Reynolds Number based on free-stream conditions 
r = recovery factor 
St = Stanton Number [Eq. (3)] 
T,)’ = local stagnation temperature 
7; = stagnation temperature as measured by a stagnation 
temperature probe 
T = local static temperature 
7. = equilibrium wall temperature for zero heat transfer 
u = velocity 
u, = friction velocity V tw/pw 
u+ = velocity parameter u/u; 
y = distance perpendicular to wall 
y* = wall distance parameter yu,/v 
y = ratio of specific heats 
= total boundary-layer thickness 
6* = boundary-layer displacement thickness 
F) 
AE er ee Jo 
J0 Polk « : 
7) 
0 = momentum thickness hes [ ei: | dy 
J0 Poll Uu. 
m = viscosity 
v = kinematic viscosity 
p = density 
T = shear stress 
Subscripts: 
LL = edge of laminar sublayer 
w = values based on wall conditions 
p = temperature probe 
0 = values based on momentum thickness 


= values based on free-stream conditions outside bound 


ary layer 


EXPERIMENTAL EQUIPMENT AND TECHNIQUES 


Wind Tunnel 

The boundary-layer surveys have been conducted in 
the NOL 12 by 12 cm. Hypersonic Wind Tunnel No. 
1% 1 11 This tunnel operates continuously in the 
Mach Number range 5 to 10 free of air condensation 
effects. Supply temperatures from 300 to S800°K. and 
supply pressures from | to 45 atm are available. A two- 
dimensional water-cooled wedge nozzle expands the air 
(see Fig. 1) and can be adjusted for any Mach Number 
simply by changing the throat opening area. 

To minimize the probe positioning errors and to 
facilitate accurate determination of the profile shapes, 
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it was decided to survey the relatively thick (approxi 
mately 25 mm.) turbulent boundary layer on one of th 
nozzle walls. Such a boundary layer is subjected to a 
slight free-stream pressure gradient in the neighborhood 
of the survey station because the flow is radial. (The 
corresponding Mach Number rise is about 3 per cent 
per tunnel caliber at the survey plane for a Mach Num 
ber of 5.0 and decreases with increasing Mach Number 

However, it is felt that the effect on the profiles is small. 
The wall temperature, on the other hand, is maintained 
practically constant near room temperature over the 
entire length of the nozzle, except at the throat where 
the surface temperature approaches the recovery tem 
perature. In all of the tests the transition from laminar 
to turbulent boundary layer occurred slightly down 
stream of the nozzle throat.* The exact effect of the 
history on the local boundary-layer characteristics is, of 
course, not fully known and the results must be in 
terpreted with this fact in mind. 


Boundary-Layer Surveys 


The boundary-layer surveys were made at the center 
of the nozzle wall, approximately 4 in. upstream of the 
nozzle exit plane (see Fig. 1). Measurements at this 
position are unaffected by the junction of the nozzle end 
and the first diffuser plate.!* Reynolds Numbers and 
the rate of heat transfer to the wall are controlled by the 
supply pressure and the supply temperature. In all 
tests the Reynolds Numbers based on boundary-layer 
thickness are roughly two orders of magnitude greater 
than needed to make slip flow effects at the surface 
negligible.'* Since the wall temperature is always 
maintained near room temperature, the lowest rate of 
heat transfer to the water-cooled nozzle that could be 
realized corresponds to a supply temperature which is 
just high enough to avoid air condensation in the test 
section. 

For each survey, pitot and static pressures, stagna 
tion and wall temperatures, and the wall temperature 
gradient perpendicular to the nozzle surface are re 
corded. 
diameter orifice in the wall and by a static probe in the 


Static pressures are measured by a 0.64 mm. 


free stream just outside the edge of the boundary layer. 
The static probe is an 8° cone cylinder with orifices 
located 17 diameters aft of the shoulder (see Fig. 2). It 
is mounted from the side wall, with its axis parallel to 
the flow. Since agreement between wall- and free- 
stream static pressures is within | per cent, a constant 
static pressure has been assumed to exist through the 
boundary layer. The measurements are made with 
silicone oil manometers of +2 microns measuring ac- 
curacy. 


* Surface probe tests made subsequent to the results reported 
in reference 13 indicated that the boundary-layer profile r 
ported in this reference was close to the transition region. It was 
definitely established, from the results of the surface probe test 
that all boundary-layer profiles presented herein were measured 


in the fully turbulent region. 
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The pitot pressure is surveyed from wall to free stream 
with a flattened hypodermic tube of 0.25 mm. opening 
(Fig. 2). The opening is large enough to avoid errors 
due to slip flow effects on the pressure measurements in 
the region of low Reynolds Number in the boundary 
layer near the wall.''| The probe is mounted in the 
microtraverse mechanism which also accommodates 
the connection to the pressure gage. The position of the 
probe in reference to the wall, can be read to +0.025 
mm. The zero position at the wall is defined as the 
position of electrical contact between probe and wall. 
Impact pressures above 20 mm. Hg are measured with 
a precision mercury manometer of +0.1 mm. measuring 
accuracy. For lower pressures the silicone oil manome- 
ter is used. 

A stagnation temperature probe with a single plat- 
inum-coated silica shield is used for the temperature 
surveys of the boundary layer (Fig. 2). For the 
measurements close to the wall, a flattened probe with 
a half opening of 0.48 mm. is used. For large Reynolds 
Numbers the temperature recovery factor of the probes 
is 0.998. The temperature recovery factor variation 
with flow parameters for each probe is described by a 
single calibration curve by relating the calibration data 
to the flow conditions inside the probe. Reference to 
this curve makes it possible to determine accurately 
total and static temperatures throughout the boundary 
layer from measured temperatures and pressures. The 
e.m.f. output of the temperature probe is recorded on 
Brown temperature recorders. The temperature meas- 
urements are accurate to +0.2°C. 
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Heat-Transfer Measurements 


Local values of the heat transfer to the nozzle wall 


and nozzle surface temperatures at the boundary-layer 


survey station are obtained from temperature measure 
ments in the nozzle wall. These measurements are 


nade with four thermocouples imbedded in the nozzle 


wall at various distances from the surface. Depending 
upon operational conditions of the tunnel and the rate 
of coolant flow to the nozzle, equilibrium readings are 
reached after 10 to 20 minutes of tunnel operation 
They show that the temperature drops linearly from the 
nozzle surface. Since previous investigations indicate 
that longitudinal and lateral temperature gradients at 
the nozzle surface can be neglected,'* the local rate of 
heat transfer can be calculated quite accurately. The 
wall temperatures are measured to within +0.01°C 
with a Leeds and Northrop K-2 type potentiometer. 


DATA REDUCTION 


The experimental data needed to evaluate the bound 


ary-layer profiles are: 


(1) Pitot pressure, po’, vs. distance from wall. 

(2) Static pressure, p, taken constant. 

(3) Measured total temperature, 7';, vs. distance 
from wall. 

(4) Variation of probe temperature recovery factor 


with flow parameters. 


Since py’ and 7’; are not always measured at identical 
positions from the wall, both sets of data are taken at 
sufficiently close intervals so that the temperatures cor 
responding to impact pressure readings at a particular 
position can be interpolated from a curve faired through 
the temperature data. 

The Rayleigh formula is used to compute the Mach 
Number profile from py)’ and p.'® All ratios of flow 
parameters needed in the further evaluation of data are 
then read from Burcher’s flow tables.'” 

To evaluate the stagnation temperature, 7)’, from 
the measured 7’; values requires reference to the calibra 
tion curve of the probe.” This curve gives the re 
covery factor of the probe r, = (7°; T)/(To’ — T) as 
a single function of Nu(k,/k,), where Nu is the Nusselt 
Number based on the flow conditions inside the probe, 
and k, and k, are the thermal conductivities of air and 
thermocouple wire, respectively. The function Nu 
(k,/k,) can also be expressed in terms of the measured 
pitot pressure and the measured total temperature. 
With the recovery factor read from the calibration 
curve, the 7)’ value can be calculated. 


y-1 


—— M?+1 


Static temperatures and velocities are then computed 
using the standard relations between Mach Number, 
static and total temperatures, and velocity." 
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TURBULENT BOUNDARY 


To calculate the factor u, for the u*, y* representa- 


tion'’ of the boundary layer velocity profiles, the wall 
shear stress 7, 1s calculated either from the slope of the 


velocity profile in the laminar sublayer 
Tw = ph (du dy 

or from the temperature measurements in the nozzle 

wall, using the following equations: 


; (AT’/Ay)k 
on = — eon (3) 
(7, — 7 ce Cpa Pp 


and 

to = Mak “6.%.* (4) 
Values of A7’/Ay are obtained from the measured wall 
temperature gradient and values of 7, were calculated 
with Pr = 0.72. 


For the cases where heat-transfer measurements were 


usmgr = Pr 


available (Table 1) the latter procedure is used to de 
termine cy. The diiferences in values of c,; obtained by 


the two procedures are in no case larger than 5 per cent. 
DISCUSSION OF RESULTS 


Boundary-Layer Profiles 


Naturally turbulent boundary layers were surveyed 
at free-stream Mach Numbers of 5.0, 6.8, and 7.7. A 
summary of the complete data is given in Table 1 and 
typical plots of the basic data are shown in Figs. 3 
through 10. 

Four surveys were inade at a Mach Number of 5.0 at 
Reynolds Numbers Re, of 5,350, 6,480, 7,950, and 7,370. 
The (YT, — Te)/Te: 


varied from 0 to 2.36 by changing the stagnation tem- 


heat-transfer parameter, was 
The surveys at lJ = 6.8 
12,640, S,400, 


perature of the air stream. 
were obtained for Re, values of 8,550, 


Pal 


Tabulation of Operating Conditions 
t 5 


LAYERS IN 


3LI 
ind Derived Parameters for all Measured Boundary-Layer Profilest 


RSONIC FLOW o 


HYPE 


and 7,960; the heat-transfer parameter was varied 
between 3.05 and 4.63. 


at M = 


5.65. 


A single survey has been made 


- T,)/T 


7.67 with Re, = 8,130 and (7, 
For most of the curves, a dimensional distance has 
been selected for the abscissa in order to show more 
clearly the physical differences between the profiles at 
different heat-transfer rates. For the same reason some 
of the curves and data points have been omitted from 
Figs. 3 through 6 and Figs. S and 9. Shape differences 
are pronounced only in the temperature profiles of the 
boundary layer, Figs. 4, 5, and 8. The static tempera- 
ture profiles obtained with heat transfer have a tem 
perature maximum close to the wall. The tempera 


tures at the maximum or those still closer to the 
wall could not be evaluated from measured air tem- 
peratures because of the physical size of the tempera- 
ture probes. In the cases where heat-transfer data have 
been measured in the nozzle wall (see Table 1), the slope 
of the temperature curve immediately at the wall has 
been deduced from the temperature gradient in the 


nozzle wall 
AT AT x 
k 


k . 
Ay/s Ay J ais 


(Values of the thermal conductivity of air have been 
taken from Eckert.” 
data and the computed slope have then been joined to 


The curve through the measured 


give a smooth shape to the static temperature curves. 
Since only the square root of the temperature enters 
into the computation of the velocity, the velocity profile 
close to the wall was found to be insensitive to errors in 
7’ due to 
measured air temperature and the wall temperature.) 


incorrect interpolation between the last 
Comparison of the temperature curves with theory may 
be utilized to evaluate the turbulent Prandtl Number. 


























| | | «| .. | neg 3° e iz | Bescispgt oo. | “ i 

| Meo P,, Atm | T,, K , ae a a. Re, 108 cgxi04 ce i. | & $1 uy = yy, 

a Cry — a | 3 

| 4.93 3.08 325 | 5.42 1.63 | 7.12 | 0.624 | 5350 | 10.9 0.3686 ty) | 0.655 | 0.0375 11.97 

5.01 | 5.01 399 | 4,29 1.65 | 6.13 | 0.708 | 6480 | 10.9 | 10.6 0.381 1,23 0.580 | 0.0272 12.44 | 

5.03 | 7.85 513 3.49 1.56 | 5.82 | 0.815 | 7950 %-43 | 9.69 | 0,341 2.07 0.570 | 0.0260 14.90 

|5.06 | 8.58 | 562 3.27 1.52 | 5.72 | 0.871 | 7370 | 9-18 = 0.329 2.36 0.575 | 0.0238 | 14.98 | 
+——_—— —— . ‘ 

6.83 | 15.3 468 6.34 2.50 | 10.91 0.784 | 8550 6.83 | 6.66 | 0.251 3.05 0.600 | 0.0440 12.90 

6.83 | 28.5 586 5.24 2.94] 9.41 | 0,903 | 12640 5.93 | 5.82 | 0.2335 4.16 0.565 | 0.0322 14,30 

6.78 | 21.1 586 5.22 2.24 | 10.34 | 0.825 | 8400 6.66 | 6.65 | 0.244 4.05 0.588 | 0.06375 14.03 

6.78 21.4 639 4.64 1.81 9.29 | 0.882 | 7960| 6.94 | 6.70 | 0.251 4.63 0.583 0.0410 14.90 

7.67 24.2 645 5.94 2.00 | 13.12 1.052 | 8130/| 5.98 -- 0.2174 5.66 0.565 | 0,049 12.95 


















































Note: Unless specified otherwise, all fluid properties involved in the Reynolds Numbers and skin-friction coefficients are evaluated 
at the edge of the boundary layer. 

+ Correction for M. = 7.67: uzt/u~ = 0.600, 62/5 = 0.058, u 

* Obtained from the velocity profile slope at the wall. 

** Obtained from heat-transfer measurements in tunnel wall. 

*** cr; obtained from Karman-Schoenherr equation for zero heat transfer 


L 


r= yp* = 14.2. 
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] Peo oe : “i Xi from the heat-transfer data using the Reynolds analogy 
1.0 on ——f TO eS — (see Eq. +). Agreement between the values of shear 
| | | . ° 
9 | ae a | | | stress calculated by these two methods is of the order oi 
. eat | | | } | 2to5percent. This result extends the range of applica 
8 ——— TT TT bility of Reynolds analogy to hypersonic Mach Num 
7 5 | | _| bers. Rubesin’s*! modified Reynolds analogy gave 
i. a | shear stress values within 1.5 per cent of those calcu 
u ‘6 a - e 4 63| | | lated by the Reynolds analogy."® 
es Pg T ' 
- ———_ «3.00 | Laminar Sublayer 
4 én — oe am weet ess ie 
| | From the results shown in Figs. 7 and 10 it is appar 
4 (_ PROFILE NEAR WALL __| +—4 ; : 
‘ Of Ss teete | | ent that the experimental u*, y+ velocity profile repre 
2 = be df sentation only roughly specifies the extent of the 
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stream Mach Number of 6.8 and two rates of heat transfer 


Preliminary evaluations of this quantity made for the 
zero heat-transfer case of \J = 5.0 using the equations 
given by Rubesin*! yielded a value for the turbulent 
Prandtl Number of the order of 0.9. 

The velocity profiles at any one Mach Number, Figs. 
6 and 9, are similar in the turbulent outer part. How- 
ever, the slope in the laminar sublayer increases with 
increasing heat transfer. 

The w+, y* representation of the velocity profiles, 
Figs. 7 and 10, show a pronounced upward shift of the 
turbulent portion of the curves with increasing heat 
transfer. A similar shift has been measured by 
ler.°2. A Reynolds Number effect is apparently super- 
imposed. This is indicated by the position of the three 
upper curves of Fig. 7 relative to each other. Accord 
ing to the respective values of the heat-transfer parame 
ters for these surveys, a different spacing in the turbu 
The concave curvature 


Deiss 


lent portion would be expected. 
of the turbulent portion of the curves is found in most 
Also shown in Figs. 7 and 10 are the 
laminar sublayer and the 
+) for the turbu- 


flat plate surveys. 
curves u+ = yt for the 

logarithmic law (u* = 5.5 + 
lent portion of the boundary layer for incompressible 
Only one survey is available 


5.75 logioy 


flow at zero heat transfer. 
for the case of zero heat transfer, and this is at a Mach 
Number of 5.0. value of the 
edge of the laminar sublayer is about 11.9 which is close 
This quantity is sig 


For this case the w+ = yt 


to the incompressible flow value. 
nificant because many of the turbulent boundary-layer 
theories utilize this value and assume it to be constant 
with Mach Number and heat transfer. 
sults indicate this to be true only for the case of zero 
heat transfer, and as the heat transfer is increased the 
value at the edge of the laminar sublayer in- 


The present re- 


ut = yt 
creases. 
The wall shear stress values needed to compute the 
friction velocity terms of the u*, y* coordinates were 
obtained either from the slope of the velocity profiles or 


closely fitted with a straight line (Fig. 11), the slope of 

which is the value of the exponent in the power profile 

representation of turbulent boundary-layer profiles, 
u/U. = (y/5)'/" (6 


Approximate values of 1 of 7, 6, and 5.5 have been 


determined for \/J = 5.0, 6.8, and 7.67, respectively. In 
the laminar caller, a straight line of unit slope may 
be drawn through the data. The point of intersection 
of these two straight lines is designated as the edge of 
the laminar sublayer. Fig. 
of this procedure applied to the data for a typical pro 
file. The intersection point so determined is fairly in 
sensitive to the slope of the line faired through the ex 
perimental data in the turbulent portion of the velocity 
therefore be determined with a fair 
It should be noted, however, that 
the thickness of the 
above that the 


11 demonstrates the results 


profile and may 
degree of precision. 
the procedure for determining 
laminar sublayer described assumes 
buifer-layer region is nonexistent. 

This method for determining the edge of the laminar 
sublayer has been applied to all the profiles included in 
The results are tabulated in Table | 
In this figure, the 
with heat 


the present paper. 
and are shown graphically in Fig. 12. 
variations of the experimental value of u,/u 
transfer parameter (7°, — 7;,)/7.. are compared with 
values of these parameters calculated from equations 
values of Re, 
analysis embodies Re; Re, 
At first glance it appears that the experi 


at appropriate 
rather than 


given by Donaldson?* 


(This as a 
parameter. ) 
mental values of uz,/u. are 
either heat transfer or Mach Number. 
examination of the data indicates that, qualitatively, 
both the trends and the values predicted by Donald- 
son's theory are found experimentally. The 
12 indicate 


relatively insensitive to 


However, closer 


results 
that 
is considered) the 


shown in the lower portion of Fig. 
(when Reynolds Number variation 
diminution of 6,/6 with increasing 
T.)/T. predicted by the equations of reference 
observed the 


values of (7, 
23 was 
seeiitel 


experimentally. However, 
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Static-temperature variation across boundary layer for a 


free-stream Mach Number of 5.0 and three rates of heat transfer. 




















25 — 
| fr 
| | // 
Te-"e . 2.36—_| i r 
. 2.07 Dag 5 A 
20} — oe ed 
u*= Uu * 44 rit genes 
> 
[Tw LP” i 
Pw [hf 
f A 
15 |_| ff” | 
Y us 5.5+*5.75 LOGY 
@ i 
| P| | 
| | 
M MEASURED, T INTERPOLATED 
— a 
| | 
50 100 500 1000 
Y*ey [Xu / 
Pw Vw 
Fic. 10. Semilogarithmic representation of boundary-layer 


profiles for a free-stream Mach Number of 5.0 and four rates of 


heat transfer 














8 JOURNAL OF 
1.00 

080} 

0.60 | 






































u 
1 4 | | 
co; it | | 
C20 1 a aa t ———— J 
O10 & = ; a oe oe 
0.0l 0.02 0.04 O10 020 040 LO 
y/§ 
Fic. 11. Log-log plot of velocity profile for a free-stream Mach 
Number of 6.8 
sai DONALDSON |e 
M* 68 THEORY nen 
[= BS Me =6.8 
ie 
u, ie, 4 Mc =7.7 
a 0.60 ren: 
Ue — ° * — i 
M=50  —~—2 . ‘ 
age 
0.50 , ; 
0.05 - 7 
ae ” , 
7 —! a 
Sawa —— s 
% = ° ee ~~! 
- a ces 
0.01 1 1 
0 2 4 6 
Te-Tw 





HEAT TRANSFER PARAMETER - 


Fic. 12. Variation of velocity and thickness ratios of laminar 
sublayer with heat-transfer parameters. 


values are somewhat lower than the experimental 
points. It can also be seen in this figure that the 
laminar sublayer thickness increases by roughly a factor 


of two as the Mach Number increases from 5 to 6.8. 


Skin Friction Results 

In the figures which follow, the experimental values of 
cy were determined by the velocity-profile slope tech- 
nique described in a previous section, and the values of 
cy, used to form the ratio c,/c;; were calculated from the 
Karman-Schoenherr equation for the same Re, value 
and zero heat transfer: 

0.02932 


= (7) 


Logi (2Rep)[1/2 login (2Reg) + 0.4545 


Cri 


The Reynolds Number based on momentum thickness 
was used in preference to the customary Reynolds Num- 
ber based on distance from the leading edge.' A 
Reynolds Number based on a boundary-layer parameter 
was selected because of the arbitrariness inherent in 
calculating an effective leading edge for boundary-layer 
measurements made on a wind tunnel wall. 

The variation of c,/c,;; with heat-transfer parameter 
(T. — T,,)/T.. is shown in Fig. 13 for both the 17 = 5.0 
and lJ = 6.8 cases. The data indicate that increasing 
values of heat-transfer parameter have little effect on 
the skin friction ratio. On the other hand, for the range 


THE AEBRONAUTICAL SCIENCES 


JANUARY, 195: 


of heat-transfer conditions investigated, Van Driesi 
and Rubesin*! predict an increase in the skin friction 
ratio of about 10 percent. It should be noted, however, 
that it is implied in these analyses that the edge of the 
laminar sublayer occurs at fixed values of w* and y 
and the slope of the turbulent portion of the u* 
curve is constant, regardless of the Mach Number o1 
heat-transfer conditions. As pointed out previously the 
results presented in Figs. 7 and 10 and Table | show 
that these assumptions are not generally valid. Pre 
liminary calculations have indicated that if the upward 


¢ 


displacement and changing slope of the u*, v* curves 
with increasing values of (7), — 7),)/7'. are accounted 
for, the value of c,/c;; decreases as the values of (7, 
7) /T.. imerease. 

The variation of c,/c;; with Mach Number is shown in 
Fig. 14. 
tion measurements of Coles.*. Inasmuch as the previous 


Also shown on this figure are direct skin fric 


figure indicated a small variation of c,/c,; with heat 
transfer parameter, it was felt that a curve connecting 
the insulated plate data of Coles” with the present data 
was ‘ustified. The dispersion of Coles’ data is a 
Reynolds Number effect and not experimental scatter. 
A similar dependence on Reynolds Number is apparent 
in Coles’ data if the Reynolds Number is based on lead 
ing edge distance. The position of the dashed curve, 
identified as Re, = Constant = 7,000, was fixed by plot 
ting the variation of c,/c,;; against Re, for each Mach 
Number picking oif the values of c,/c;; at a value of Re,= 
7,000. This value of Re, was selected because it repre- 
sented a Reynolds Number for which the greatest over- 
lap of data existed. It is apparent that the c,/c,; varia 
tion with Mach Number is a smooth continuous curve 
for a single Reynolds Number. 

From the results shown in Fig. 14 it may therefore be 
concluded that the present data are in good accord with 
Coles’ data.’ The skin friction results presented are 
consistent with the trends predicted theoretically; 
the numerical values are in good agreement with those 


predicted by Wilson’. 
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CONCLUDING REMARKS 


Detailed explorations of turbulent boundary-layer 
velocity and temperature profiles have been made on a 
nozzle wall of the NOL 12- by 12-cm. Hypersonic Wind 
Tunnel at Mach Numbers of 5.0, 6.8, and 7.7 for vary 
ing rates of surface heat transfer. From these meas 
urements it has been found that turbulent boundary 
layer profiles in hypersonic flow qualitatively resemble 
turbulent boundary-layer profiles in incompressible 
flow in many details. In particular, all of the velocity 
profiles measured can be fitted with a power law in the 
outer turbulent portion. As in the case of incompressi 
ble flow, the turbulent portion of the profiles differs in 
shape from the logarithmic-law velocity profile. In 
this region the discrepancy between the experimental 
data and logarithmic velocity profile law becomes 
greater with increasing heat transfer. 

It was found that the velocity proSiles in the laminar 
sublayer are linear. Furthermore, the ratio of the 
laminar sublayer thickness to the total boundary-layer 
thickness decreases slightly with increasing heat-trans- 
fer rate but increases considerably with increasing 
Mach Number at common heat-transfer rates. The 
u* = y* value at the edge of the laminar sublayer for 
the zero heat-transfer case is close to the incompressible 
flow value. However, this w+ = y* value is found to 
increase with increasing heat transfer for any given 
Mach Number. 

It is demonstrated that the Reynolds analogy is ap 
plicable for turbulent boundary layers in hypersonic 
flow. The values of skin friction coefficient calculated 
from the heat-transfer measurements by using the 
Reynolds analogy agree to within 5 per cent with those 
skin friction values deduced from the velocity profile 
slope in the laminar sublayer. These values are found 
to be only slightly affected by heat-transfer rate and 
are in accord with the direct skin friction measure 
ments of Coles when the results are based on a single 
value of Reynolds Number. 

The present results extend the range of available skin 
friction data toa Mach Number of 7.7. Finally, the de- 
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tailed measurements of velocity and temperature dis 
tributions across turbulent boundary layers in hyper 
sonic flow presented herein enlarge the fund of available 
data for turbulent boundary layers in compressible 


flows. 
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A Study of Surge in Compressors and 
Jet Engines 


CARL E. PEARSON* 
Harvard Unwersity 


SUMMARY 


Some experimental results concerning compressor surge are 


presented. A theoretical surge model suitable for either a com- 
pressor or a jet engine is constructed, and the equations governing 
its behavior are analyzed. Certain surge criteria are deduced, 
and their implications for both steady-state and transient oper- 


ation discussed. 
(1) INTRODUCTION 


Pree YRS AND JET engines are occasionally sub- 
ject to a form of instability called ‘‘surging,”’ 
which usually occurs near the peak of the pressure-flow 
characteristic and which seems to involve an overall 
oscillation of flow. Some general surveys and some 
analyses concerning this phenomenon will be found 
in references | to 6. 

The present analysis deals with a theoretical model 
of a jet engine (or compressor, as a special case) con- 
sisting of a compressor, a “‘burning’’ chamber in which 
heat but not mass is added, and an exit nozzle which is 
assumed choked so that the turbine has no aerody- 
namic coupling to the compressor. Shaft coupling be- 
tween turbine and compressor is assumed to have but 
a secondary effect on the incidence of surge. A num- 
ber of idealizations (such as the assumption that 
burner temperature varies linearly from inlet to exit) 
are made in order to render tractable the manipulation 
of the mass and energy conservation equations as 
applied to each element of the model. A certain ex- 
periment indicates that the character of the airflow in 
the compressor under transient conditions depends 
only on the instantaneous flow rate and compressor 
r.p.m.; this generalization allows a type of entropy con- 
servation theorem to be written for the compressor 
section. 

The resulting set of ordinary differential equations 
must be obeyed by the airflow at all times; it is now 
necessary to examine these equations to see if solutions 
exist which oscillate with increasing amplitude. Be- 
cause of the nonlinearity of the equations, a direct in- 
stability analysis is not feasible; therefore it is assumed 
with reasonable justification that if a solution becomes 
unstable then so does any perturbation in this solution. 
The instability of the set of linear equations resulting 
from perturbation of the original equations is examined, 
and three analytic criteria obtained. 

Received February 17, 1954. 
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Within the limitations of the model, the following 
facts concerning surging may be deduced from these 
criteria: 

(1) A jet engine surges as a unit—.e., all elements 
of the engine experience the same oscillation in airflow. 

(2) A jet engine surges at approximately the same 
operating point as its compressor when tested alone. 

(3) Surge tends to occur when the slope of the 
throttle characteristic is greater than that of the over 
all compressor characteristic. 

(4) The surge point of a compressor is influenced to a 
small extent by the volume of the surge chamber, in 
such a manner that although a compressor will surge 
at a certain point when tested with a large surge cham- 
ber, it may not surge in a jet engine at the same oper- 
ating point. 

(5) Any alteration in the engine which increases the 
losses due to boundary-layer separation will tend to 
make surge occur earlier. Such losses might be ac- 
centuated by an increase in the number of diffuser 
vanes (except in those cases where the increased num 
ber of vanes inhibits separation) or by operation at 
sufficiently low Reynolds Number (as at high altitude) 
that certain boundary layers lose their turbulent char- 
acter and so separate earlier. 

(6) A large rate of increase of burner or surge chamber 
pressure may tend to induce surge if it occurs in the 
neighborhood of a critical operating point. Suchincreases 
would occur in the jet engine case if the rate of fuel 
injection were suddenly increased (e.g., by opening of a 
secondary manifold) and in the compressor case by 
running the compressor into surge along a constant 
speed line rather than along a constant throttle line. 


(II) SomE EXPERIMENTAL RESULTS 


The present work was carried out on a centrifugal 
compressor consisting of a separate axial flow inducer 
and centrifugal impeller, the spacing between the two 
units being variable. A description of this apparatus, 
and a survey of its surge and stall-propagation char- 
acteristics, will be found in reference 7. 

It was considered that the behavior of the inducer 
was representative of a single stage of a multistage axial 
compressor, the impeller here providing an effect equiva- 
lent to the subsequent stages in an axial machine. 

Beyond the results of reference 7, one rather inform- 
experiment involving transient operation of 
One hot-wire was 


ative 
the compressor may be described. 
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SURGE IN COMPRESSORS AND JET 


placed just ahead of the inducer, and the voltage across 
this wire was amplified by a direct-coupled amplifier 
and presented on an oscilloscope. Holding the r.p.m. 
constant, the throttle was set at various positions and 
the oscilloscope deflection noted. A second hot-wire 
was placed between the inducer and impeller; the volt- 
age across it was diiierentiated, amplified, and pre- 
sented on a second oscilloscope screen. By means of a 
mirror system, both oscilloscope traces were super- 
posed and arranged for photography, using a high- 
speed continuous-strip camera. The compressor was 
then set at a certain operating point for which the 
throttle was wide open; the throttle was then suddenly 
closed and the responses of the two hot-wires were re- 
corded. The film speed was such that the (instan- 
taneous) blade wakes could be examined with ease. 
The compressor speed did not change appreciably during 
the '/1o sec. it took for the flow to alter from wide open 
to deep surge. It is unfortunately not feasible to re- 
produce the film strips resulting from such experiments, 
but the main result may be concisely stated as being 
that the instantaneous character of the air flow in the 
inducer blading during transient operation was ex- 
actly the same as that corresponding to a steady mass 
flow equal to the instantaneous value of that flow. 
This result indicates that the transient operation of 
compressors may be analyzed on the basis of their 
steady-state behavior; some further evidence toward 
this effect is to be found in the rapidity with which the 
character of a boundary layer accommodates itself to 
a change in the angle of incidence, as demonstrated in 


the phenomenon of stall propagation (reference 7 
Although this conclusion concerning the identity of 
steady-state and transient airflow character can be only 
approximately true for a compressor involving a large 
number of blade rows, its use should not alter the gen- 
eral character of the results to be deduced from a study 
of a surge model based on its adoption. 


(III) SuRGE MopEL 


The geometry of a jet engine will be modeled by 
means of an enclosure, representing the effective burner 
volume, connected to inlet and exit channels. The inlet 
channel contains the compressor, which is assumed to 
draw air from the atmosphere. The effect of com- 
bustion in the burners will be obtained by the addition 
of heat from an exterior source. The exit channel will 
be considered choked. In the case of a compressor 
test, the burner volume becomes the volume of the 
surge chamber, and of course there is no heat addition 
whether the exit channel is choked or not has little 
effect on the general character of the analysis). 

The arrangement is shown in Fig. 1, and the following 
notation will be used: 


A = area of channel, ft.? 

C = constant 

C, = (constant) specific heat of air at constant pressure, ft.- 
Ib./lb. °F 
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Fic. 1. Surge model 

C, = (constant) specific heat of air at constant volume, ft.- 

Ib. /Ib. °F 

dx = element of channel length, ft 

dV = element of volume, cu-ft 

g = acceleration of gravity, ft./sec.? 

H = rate of heat addition, ft-lb. /sec 

K = constant 

P = pressure, lb./ft.? 

P, = atmospheric pressure 

P, = pressure in enclosure 

Q = inlet channel constant, ft. sec.?/Ib 

R = gasconstant = P/rT 

r = density, lb./ft.* 

Yav. = average density in burners 

r, = density of surrounding atmosphere 

r, = density at point 1 

S = specific entropy, ft.lb./lb. °F. Abs 

TY = temperature, °Abs. Fahr 

7, = temperature of ambient atmosphere 

T, = temperature at point ‘1,’ Fig. 1 

7, = temperature at point ‘‘2,’’ Fig. 1 

t = time, sec 

V. = volume of burner enclosure, ft.* 

V; = volume of region 7 of Fig. 1, ft.* 

W = flow rate, Ilb./sec 

W, = value of W in compressor channel 

W. = value of W in exit channel 

@ = compressor channel specific entropy change, ft.lb./°I 


Ib 


= work per lb. of air done by compressor, ft.lb./Ib 


Ww 


Consider first the compressor channel. The require- 
ment of conservation of energy applied to the region 
t of Fig. 1 bounded by the dotted lines ‘‘0” and “L”’ 


yields the equation 


WiC,T, + wW: — W,C,T; = 


of es W,? ; 
C,T + rd\ 
ol eo (1) 2¢er°A? 


where the kinetic energies of the air at surfaces ‘‘0’”’ and 
“1” are assumed negligible. Since not much air can 
accumulate in the compressor blading, W; is approxi- 
mately constant over the region 7 and hence the last 


integral may be rewritten as 


Q W? ow. 
—_ rA dx = Wi — Q, 


Ot 2er? a* ot 
, ‘ dx 
where the channel constant Q, is given by ‘i —”; 
v(t) 2s 


may be considered as a constant because of the 
fact that the variation in density is negligible com- 
pared with the variation in W,. Using this result and 


also the gas equation, 
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Finally to avoid the awkward term | P dV, let us 


approximate it by 


‘  & 
PdV =~ (Pp + P,) 


« ~ 


where A is a constant of order unity. This approxi- 
mation will turn out to be not too important; in any 
case, it should not affect general surge tendencies. 
The energy equation for the compressor channel may 
therefore be written: 

KC,V; Pe , OW, 


set QO, (1) 
2R OW, Of 


CT. C,T) + w= § 
where P, = OP,/dt. 
We require also the energy equation for the enclosure : 


a) 
WiC,Ti + H = WiC,T: + > CrT dV 


e (e) 


. 
where | indicates integration over the burner vol- 


ume, and where again kinetic energies may be neglected 
in comparison with the remaining terms. Using the 
gas equation, 
- es CVes 
WiC,7, + HW = W2C,T2 + Pp. 2) 
R 
where we have made use of the fact that (even in tran- 
sient operation) the pressure has practically the same in- 
stantaneous value throughout the burners. 

The compressor channel has another property of 1m- 
portance in addition to that involving energy, viz., 
the losses introduced by the blading, channel walls, 
etc. These losses are best handled by making use of 
the entropy function S = C, In 7 — Rin P. It is 
shown in Appendix I that an entropy conservation 
theorem may be written for the compressor channel; 


this equation takes the form: 
W W,| Cyn — RI = + 
1p = 1 j p inl r, n P. 
nF .. as , 
| (C,InT — Rin P)rdl 
Of J (i) 


The order of magnitude of each term in this equation 
is easily calculated, and the result is that the last term 
may be neglected without appreciable errors. Hence 


2 


je P, 
r, + RI1n P, = 0 (; 


@+ C,1n 


Mass conservation in the burners requires that 


, , a) ; 
W, — Ws = * f rd\ 
of (e) 


Now the integrand of this last expression is going to 
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SCrENCES— 
depend on the manner in which heat is added. The 
simplest reasonable assumption which can be mad 
here is that the average temperature through the 


burners 1s 


and therefore 
ra 2 
R 74+ 7: 


The mass conservation equation becomes 
W, — We = Vela 5 


Finally, since the exit nozzles in a jet engine are 
normally choked, 
P 
We2=C 6 
VT> 


where C is a constant of proportionality. Strictly 
speaking, this last equation should take into account 
time-dependent effects, but in the simpler case of a 
nonchoked nozzle it may be shown that these terms 
are comparatively small; hence we will omit them here 
also. 

The six numbered equations involve the six unknowns 


Wig Ws, 2 yy De Pe Kan 


and we therefore have a complete set of equations. 
These equations are nonlinear and therefore difficult 
to handle; although it is true that it is possible in prin 
ciple to combine all these equations into one by suc 
cessive elimination, the order and complexity of the 
resultant equation are such as to make it impractical 
to use graphical phase-space analysis (this being the 
usual technique of handling nonlinear stability prob 
lems). Further, such an approach, even if feasible, 
would not yield analytic surge criteria. The phi 
losophy of solution of these equations is discussed in 
the next section. 
(IV) SOLUTION OF EQuaTIoNs (1)-(6) 
The equations may be listed as: 
5GY) P, OW, — 


C,T, — C,T; +o - _ 
; iis ( 2R / W, n= 


IW4iC,T,. + H — W2C,T> — (Ss ‘ P, = 0 
@ + C,1n (T,/7T1) — R1n (P,/P,) = 0 
ve ‘“ 2 
(Ti + Ts) rave — (3) P,=0 
W, — We — Veay. = 0 


WiVTe — CP, = 0 
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Consider now an engine in either steady or transient 
operation. Then at every instant, the values of IV, 
o, P., r.p.m., ete., will be determined by the dynamic 
characteristics of the engine (which involve the turbine 
properties, moment of inertia of the rotor, etc.); in 
particular, the above six equations (which deal only 
with the compressor and burners) must be satisfied at 
every instant. A strict instability analysis would re- 
quire us to add to these six equations the appropriate 
remaining equations involving the other portions of 
the engine; however, the experimental evidence is that 
the rotor r.p.m. does not alter appreciably during the 
onset of engine surge, and since this is the only way by 
which the turbine and compressor may interact for 1n- 
stability we are justified in omitting from consideration 
these other equations. Even stronger evidence for 
this reasoning is afforded by the fact that an engine 
will, in general, surge at the same point as its com 
pressor would if tested alone. Should the reader still 
have doubts regarding this approach, he can consider 
the above equations as being the only ones involved in 
those perturbations for which engine r.p.m. does not 
change—-t.e., any stability criteria which are derived on 
this basis are necessary but not sufficient. 

Consider then, say, the acceleration of an engine. 
Each quantity involved in our six equations must be 
such a function of time that the equations are satisfied. 
During the engine acceleration, surge may or may not 
occur--experimentally, the borderline is often slender. 
It seems reasonable to suppose, even if surge did occur 
during an actual acceleration, that there existed a 
surge-free set of solutions of the dynamic equations 
governing the engine, but that this solution path was 
unstable. Even if this assumption is not strictly true 
throughout the acceleration, it seems certain that a 
nonsurging solution could exist for at least a short dis- 
tance beyond the point of surge incidence, this solution 
becoming, however, unstable at this point. This is the 
assumption on which we will base our further analysis; 
the fact that the resultant surge criterion is physically 
very plausible will add confidence to the method. 
Then a perturbation approach is applicable; let us write 
each quantity as the sum of a “‘basic”’ and a perturba- 


tion portion—thus 
Wilt) = Wilt) + W(t), etce., 


where the bar over the letter indicates average value. 
Using the fact that 


du 
wee ee W,’ 
no 

: Op 
= ee W,’ 
Se. a 


where in accordance with the first paragraph the 
r.p.m. dependency is neglected in the perturbation only) 
we obtain, after linearizing, 
CT! Ow We," (*5 ) l “ 
~&y 7. " a — 
© OM | 2R / Wi 
(WiP.’ — P.W,') — Q: (OW,,/dt) = 0 


— 0 
R 
ra) (, R 
A MW a — yg as P 0 
Ol; 1; P 
1 + 7 } i. F — 0) 
R 
W,’ 1] Vi (0) 
i" 
We VT ——— 7 CP,’ =0 
yA £ 


This set of six linear equations is now easily exam- 
ined for stability, provided that the “‘basic’’ quantities, 
Il, P, ete., which are, of course, time-dependent, do 
not alter appreciably during the onset of surge—and 
this is indeed the situation experimentally observed. 
The details of the Routh-Hurwitz procedure as applied 
The fol- 


where for simplicity the bars have 


to these equations is found in Appendix IT. 
lowing conditions 
been dropped again —are found to be necessary (and in 


practice also sufficient- see Appendix) for no surge: 


Ow _. OD | I” ( =) 
—- +n i+K > (1+— )}b- 
Ol, ol, | 21, 2W,/) 
, KCV, O,R peso: OF aoe 
“ORW? | CV.P, ids ore 
(A) 
Pe 
W272C, 4 — (J, + 7 4 
2 17 
K V,C,7 
> & 
2V. W, 
Ow f Wy R 7; | Od 
" om, | We Cy 7 ' OW, 
KC,V [ WRT, | — RT; (: ; 
>. od ) 
RW? ind Ws TJ)” 
Ow 472, C (: )] W, 
aw, V7; | 2C, ] ” W. 
C ( J 4} 
i ~ W./s 
_, Oo fT> E Ce (, ay 
‘aw, (7, L° * 2C, T. 
14K 2 y (; ny 
, 7 ee T,/ § 
(C) 


> KG,V, LE  C¢ x (1 ")'] 
*ORW,2 \7,|° ° 2G, i 


Ws . G ( W, 
2 a 1+ 2 
W. ce Vs 


ve |. x j ( z) 
(7 To){ 1 + 1 + = 
a os ( T.) * T, 
1 (: 7s ) x (2 rs W x ) > 0 
V Sev: T: ee a 


It will be noticed that all these conditions are of ap- 
proximately the same form, and also are physically 
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realistic—for example, a positive 0w/OW,, a negative 
0¢/0IV;, and a positive P, all tend to induce surge— 
as they should. Consider as a special case a com- 
pressor in steady-state operation; then P, = 0, Wi = 
W., and if H = 0 then Eq. 2 requires that 7; = 7». 
Then the numerical quantities involved are usually 
such that these conditions can be (very approximately) 
replaced by the single condition 

D if Ow r Op TR 
(9) . + 


- 0 
ow ow . 


But this is precisely the condition that the slope of the 
throttle characteristic be greater than the slope of the 
compressor characteristic—for if the subscript ‘‘0” de- 
notes stagnation conditions, 


Ow Op ee: * ee 
a a i a : i = C,T 10) 
ow aw Ww we 
ae “a 
To a (C,1n Ti — Rin Po) + 7 WV 
and therefore 
iF «4 IP « | 
ee ed > 0 for stability 
dw | compressor dV | nozzle . 


This rule of thumb, condition (D), is physically very 
reasonable, and is therefore a strong indication of the 
legitimacy of the various approximations made in the 
analysis. Condition (D) is, however, very crude, and 
moreover does not exhibit the effect of transient opera- 
tion nor the effect of altering the design parameters. 
Conditions (A), (B), and (C) will be examined in the 
next section. 

Before proceeding, however, a remark or two con- 
cerning nonchoked nozzles is in order. The above an- 
alysis assumed that the nozzles were choked, and 
neglected the kinetic energy of the airin the nozzle. If, 
as is customary on a compressor stand, the nozzles are 
not choked, then a similar analysis may be carried out. 
Conditions (A), (B), and (C) are again obtained in 
similar form—the final terms are replaced by terms 
proportional to the slope of the nonchoked throttle 
characteristic, and the coefficients of the remaining 
terms are slightly altered. In this case it is easy to in- 
vestigate the effect of the kinetic energy term in the 
throttle—i.e., a second channel constant Q2. may be 
introduced—-and when typical numerical values are 
substituted in the result it turns out that this term 
is quite negligible. We were therefore justified in neg- 
lecting such effects in the choked nozzle analysis 
a result which we might also deduce from the fact that 
the Q; terms are of small order in conditions (A) and 
(C) and do not appear at all in (B); and certainly Q» 
has a much smaller effect than Q). 


(V) DIscUSSION OF RESULTS 


It is probably clearer to leave the stability conditions 
in the form of (A), (B), and (C), although the reader 
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may, if he wishes, replace 0w/OW, and 0¢/0W;, by their 
expressions in terms of total temperature T\) and 
pressure Pio as measured under steady-flow conditions 
at the compressor exit: 


Ow _ . OT 1» 

ow, "ow, 
Od = “ OT 10 RT 9 OP 19 
ow, ’odW, Py» OW, 


The following remarks may be made: 

(1) The solution of the perturbation equations is 
such that when surge commences all portions of the 
model are affected, with the same frequency but possibly 
with different phase. 

(2) Condition (D) shows that, at moderate to high 
pressure ratios (where these terms dominate), such 
things as burner heat addition and ducting volumes do 
not appreciably affect the onset of surge—i.e., the surge 
points of a jet engine and of its compressor are approxi- 
mately the same. 

(5) In general, surge tends to occur when the slope 
of the compressor characteristic exceeds that of the 
nozzle characteristic, as shown by condition (D),. 
This general statement must however be modified by 
the more exact requirements of conditions (A), (B), and 
(C); the modifications are particularly important at low 
or medium pressure ratios, such as are met during 
accelerations. Note that at high pressure ratios the 
curvature of the characteristic is usually so sharp that 
irrespective of the actual slope at which surge com- 
mences, the surge point is for practical purposes at the 
peak of the characteristic. 

(4) If a compressor stand has a large surge chamber, 
so that V, > V;, then the last term in condition (A) 
is negligible and the compressor will therefore violate 
(A) and so surge at the peak of its characteristic. For 
intermediate values of V,, the surge point will move 
toward the lower flow region; at low pressure ratios, 
the surge point of a compressor may be appreciably 
affected by the surge chamber volume. This may be 
the reason for the apparent ability of a jet engine to 
sometimes accelerate through a surge region without 
surge, as the compressor surge region may have been 
determined by use of a large surge chamber. 

(5) All the criteria show that any modification of 
compressor operation which decreases 0¢/O]V; in criti- 
cal regions tends to induce surge. Such a modification 
would be, for example, the widening of a diifusing 
throat or the insertion of extra diffusing vanes, either 
of which may increase the amount of separation losses; 
another would be the lowering of the inlet Reynolds 
Number (altitude effect) in those cases where this 
lowering destroys the turbulent nature of boundary 
layers. 

(6) At medium pressure ratios, the terms 0w/OIV; 
and 7; 0¢/0W,, no longer dominate all others; thus the 
OP,./ot term may be quite important. On a com- 
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pressor stand, there may be a difference in the surge 
characteristics depending on whether the compressor 
is run into surge along a constant throttle line or along 
a constant speed line; in the latter case, there is a 
tendency for the inertia of the air in the compressor to 
pile up air in the surge chamber and so result in a posi- 
tive OP,/dt as the throttle closes, thus inducing surge 
at an earlier point by virtue of condition (A). In an 
engine, the primary effect on OP,/Ot is as a result of 
fuel injection; if there is at any time a sharp increase 
in fuel flow rate, then OP,/Ot may have a large positive 
value and so contribute to the occurrence of surge. 

(7) The terms involving Q; and 7)/7». have only a 
minor effect on the surge point. The Q; terms are such 
that at low pressure ratios an axial compressor is less 
subject to surge than a centrifugal compressor at a com- 
The ratio 7\/7>2 is usually 
fixed by other considerations and is of little interest 
except in special cases. It may be noted that in the 
case of a centrifugal compressor, 0w/OIW;, is often ap- 
proximately zero, and the stability criteria are conse- 


parable operation point. 


quently simplified. 
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APPENDIX I 


CONSERVATION OF ENTROPY IN A COMPRESSOR 


Consider any closed surface instantaneously fixed 
in space. If the outward unit normal is m and if the 


fluid velocity is V, then Gauss’ Divergence Theorem 
applied to the surface yields for the specific entropy S: 


| Sr V-ndA = f grad:-(rSV) dV 
A Vv 


= f (r grad S:V +S grad-(rV) dV 
V 


Continuity requires 


or 


svad -(rV’) ~ 
a or 


and from the definition of material derivative 


dS oS : 
= — + grad S-| 
dt ot , 


Hence 


oa */ dS os or . 
SrV-ndA = r - 7 — § —)d\l 
JA V dt ot ot 


1.€., 
A eo - 
rS d\ (¢) 
t Jt 


dS . : 
raV = Sr V-ndA 4 
Jv dt A 


From the experimental results quoted, the left-hand 
side is a function of the instantaneous velocity alone 
and may therefore be obtained by steady-state meas- 
urements. Eq. (3) now follows—the fact that a rotor 
is contained within the region 7 does not invalidate 
Eq. (7). 

(It turns out, however, that this result is actually not 





of much present value, because in Eq. (3) the 0¢/0/ term 
was found negligible and omitted.) 


APPENDIX II 
SOLUTION OF PERTURBATION EQUATIONS 


The set of six perturbation equations given in Section 
IV must, by theory of such equation sets, have solutions 
of the form e“, or more generally of the form of a 
polynomial in ¢ multiplying e“. In any case, the 
growth or decay of any perturbation depends on the sign 
of the real part of a and this can be determined by 
examining the following determinantal equation satis- 
fied by a (where for simplicity all average value signs 
are omitted): 


Q re Ow ~ KC, V; P 0 C 0 KC,V; 0 
™~ * a, see" ° RW, * 
OP —C,T» WC, —WU oC, — R Qa 0 
oO te e 
° 0 — — 0 0 
0 = OW, T; P, 
9 
0 0 ™ ™ — — T 7. 
r r, R 1 + 
l —] 0 0 0 — V.a 
- W. : 
0 VT? 0 —( 0 


2VT» 
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which may be reduced to 


E- _. OO KC,V, | AC,V, TiR 
aQ; — 0 


— 7 P. ¥ - 
OW, ov, | RW ~ 2RW, P, 
WW, + We) 2 2% 4 or ( a) . | ew + wf in+? 
(Z 59 ») 2 a — 4 Z - ) — ( 4 
lies . *" C, WW, | C,R CP. ee 
7 W. WC 
“ite? Beye = = (7, + r | 
P. C,T2Rray. 
WiT\ Op a ar WoT> W,7,R 
+1, + T: -2 + — T2V,0 
c a, *" P, P.C, 


This determinant may be written in polynomial form: 


F303 + Fra? + Pia + Po = O ay 


where iv 
r, = 20, V2 ; 
a ee . 
pt 
2V.2T2C, dw _. Ob | _¥ ( me . ECV, I 
F, = _ 1 1 + kf 1 + — P. + 
C,R | OW, T ‘ow, | ” OV, 21V,/$ 2RW? d 
O.R ye... 3 * Ca... EK VV; C7 
ae )_WAT.R 1 Cy, + — wt. + — 74+ T8 Bde ee A fs ‘ 
C,V.P, | 2 2 (7. f 2v.N st 
BPP ow (7: E 4 Cc. (: ry | 2”, .C ( my) - 
f= = ) t : + = +t 
P, — OW, UT; 2C, Ty; W. . WW. /§ aut 
a) 1 ( 7;\? KC, V T gt 
y SS. 2 E ( te) [+4 : (: 4 yj fk 
OW, V7 2C, T2/ | i? T,/ § ise 
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Some Applications of Generalized Harmonic 
Analysis to Gust Loads on Airplanes 


HARRY PRESS* ann JOHN C. HOUBOLT? 
Langley Aeronautical Laboratory, NACA 


ABSTRACT 


A review is first made of the various aspects involved in the 
application of the techniques of generalized harmonic analysis 
to airplane dynamic behavior in rough air These aspects 
include: the concept of a power spectrum for the compact 
characterization of a random disturbance in time, the use of the 
spectrum for calculating airplane response to gust disturbances, 
and the relation between the statistics of the time history and the 
power spectrum A recently developed technique due to J 
lukey for the derivation of spectra estimates from time history 
data is also described 

The available information on the spectrum of atmospheric 
vertical velocity is summarized and a representative spectrum 
selected for use in the applications 

Some of the significant applications to gust load problems made 
to date at the NACA are then presented 
are intended to illustrate the effects of three important factors 


These applications 


involved in the airplane response to gusts, namely, the center of 
gravity position, short period damping, and wing bending 
flexibility In each case, experimental verification of the cal 
culated results is given 


INTRODUCTION 


— THEORETICAL consideration of airplane behavior 
and loads in rough air, in which the continuous na- 
ture of the rough air is explicitly considered, has until 
recently received little attention. This situation has 
existed because analytical treatments of random type 
disturbances have been too cumbersome to be practical, 
and further because information on the velocity fields 
of atmospheric turbulence was lacking, due mainly to 
the difficulty of measuring these velocities. As a conse- 
quence, almost all reductions of airplane measurements 
to obtain gust data and most analyses or airplane re 
sponse to rough air have been based upon the response 
to simple and discrete gust disturbances. This proce 
dure has also been followed in the development of de- 
sign load calculations. Although the approach has 
proved useful, it was obviously a simplification whose 
area of applicability was restricted to certain classes of 
airplanes. A need therefore existed for more refined 
and more generally applicable techniques. During the 
last few years, considerable progress has been made in a 
more complete description of the characteristics of the 
atmospheric turbulence velocity field and in the analysis 
of airplane behavior and loads in continuous rough 


Presented at the Aeroelasticity Session, Twenty-Second 
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air.'~+ This progress has been made possible by the 
development of mathematical techniques for the 
treatment of random functions by means of generalized 
harmonic analysis.° Because of the applicability of 
these techniques to gust load and to other aeronautical 
problems, it seemed appropriate to present a paper re- 
viewing the general technique and some recent gust load 
applications. 

The purpose of this paper is to briefly review several 
fundamental aspects of generalized harmonic analysis 
that are useful in aeronautical applications and to 
present some results obtained at the NACA in several 
applications of these techniques to gust load problems. 
The aspects of the theory to be considered are 

(a) The concept of an energy or power spectrum. 

(b) The use of the spectrum in relating the airplane 
response to the disturbance. 

c) The relation between the spectrum and various 
statistical characteristics of the time history. 

(d) Methods of estimating power spectra from ex- 
perimental data. 

These techniques are applied in several sets of calcu- 
lations of airplane responses to continuous rough air. 
The results obtained indicate the effects on gust loads 
of airplane center of gravity position, low damping in 
short period, and wing bending flexibility. In each 


case, experimental confirmation is also given. 


SYMPOLS 


k = statistical number of degrees of freedom, Eq. (16) and 
(17 
m = number of uniformly spaced points over the frequency 


range of interest at which power estimates are 
desired 
n = number of equally spaced observations taken over the 


record length of the time history, 7'/ Af 


Fr = period of short period oscillation, sec 
§ = stress 
t,7 = time, sec 
V = air speed, ft. per sec 
y(t) = random function of time 
yl) = derivative of y with respect to time 
y"(t) = second derivative of y with respect to time 
At = time interval, se 
An = acceleration in incremental number of g’s 
i Pa 
R(r) = autocorrelation function, lim | | y(t)y(t + 1) dt 
—aolle 
o = root mean square value, VY y2(; 
w = frequency, rad. per sec 


Q = reduced frequency, w/V, rad. per ft. 
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F(w) = Fourier transform of y(t), Eq. (2) 
#(w) = power spectral density function, Eq. (3) 


A bar above a quantity designates the time average of the 


quantity. 


GENERALIZED HARMONIC ANALYSIS 


As a preliminary to the application of the techniques 
of generalized harmonic analysis, a brief review is pre- 
sented of some of the fundamental aspects involved in 
practical applications. More detailed treatment of the 
various aspects covered can be found in references 5 to 8. 


Power Spectrum 

In the analysis of arbitrary disturbances, a frequency 
analysis by means of classical harmonic or Fourier 
techniques is frequently useful. However, for dis- 
turbances of a nonrecurring and nondissipating nature, 
these conventional techniques do not apply. On the 
other hand, for a certain class of these nonrecurring dis- 
turbances, the “‘stationary random’”’ type (that is, dis- 
turbances whose average characteristics do not change 
with time), a frequency analysis is possible by means of 
the so-called “‘power spectral density function.” To 
define the power spectral density ®(w) of a random func- 
tion of time y(¢) let 


yr(t) = y(t) -TeysT 
tr(t) = 0 elsewhere (1) 
The Fourier transform of y;(¢) always exists and is given 


by 
F(a) = | yr(te'™ dt (2) 
The power spectral density is then defined as 


l 
.|F(w)|? = lim — F(w)F(—w) (3) 
! 2a 


[To = 


P(w) = li ; 
1“ bee 271 


where the bars indicate the modulus of the complex 
quantity F(w). A point to note is that ®(w) is real and 
positive, and that unlike conventional Fourier repre- 
sentations phase is lost and only amplitude is retained. 
From Eqs. (2) and (3) 


1 ox x 
P(w) = lim 5; f f yr(byr(ee'e de dt (A) 
T—o 2r71 « —" F 


Letting « = ¢ + 7 and interchanging the order of inte- 


gration yields 


> ies l : 
_* J | ee J yr(t)yr(t + 7) dt |e” dr 
TJ) —~oLT—?o = 7 — 


(5) 
Se -; . 
= R(r)e *" dr (5a) 
ts —<« 
where R(r) is the autocorrelation function defined by 


1 ¢7 
R(r) = lim f y(t)y(t + 7) dt (5b) 


T—>« 2T. —T 
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It can be shown simply that the autocorrelation func 


tion R(r) has the following properties: 
(a) R(r) = R(—7) 


(b) R(O) = y(t) > O 
(c) R(O) > R(7r) 


From Eqs. (3) and (5a), it follows that ®(w) is always 
positive and is likewise symmetrical about the origin. 
Using these symmetry considerations and the Fourier 
integral theorem, the following reciprocal relations be 


tween ®(w) and R(r) are obtained 


of 
P(w) = } R(7r) cos wr dr | 
TS) 0 
. (6 
R(r) = | @(w) cos wr dw 
These reciprocal relations are useful in practice for ex- 
perimental determinations of spectrums as will be indi 


cated. From the foregoing it is clear that 


R(O) = y(t) = o? = | P(w) dw (7 
where a, the root mean square value (or standard devia- 
tion) of the random functions, has been introduced for 
brevity in notation. The quantity y’(¢), or 0”, the time 
mean square, provides a measure of the disturbance 
energy per unit time and has thus been characteristic- 
ally termed the power, as a carry-over from its early 
application in the fields of communications and tur 
bulence, where it often had the dimensions of power. 
¢(w) has in turn been termed the energy or power spec- 
trum. In this form, the element ®(w) dw gives the con- 
tribution to the mean square of harmonic components 


of y(t) having frequencies between w and w + dw. 


Relation Between Disturbance and Response Spectra 

A particular useful and simple relation exists for 
linear systems between the spectrum of a disturbance 
and the spectrum of the system response to the dis- 
turbance. The relation is derived, for example, in 


reference 7 and 1s 


Po(w) = &;(w)T?(w) (S 
where 
Po(w) = the output spectrum 
®;(w) = the input power spectrum 
7(w) = the amplitude of the admittance or fre- 


quency response function which is defined 
as system response to sinusoidal disturb- 


ances of various frequencies 


Eq. 8 indicates that the response at a given fre- 
quency depends only on the input and the system ad- 
mittance characteristics at that frequency which is 
readily plausible for linear systems. The application of 
this relation to the gust loads case is illustrated in Fig. 
1. The top sketch on this figure is the input spectrum 
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and in this case represents the spectrum of atmospheric 
vertical velocity. The frequency argument 2, which is 
27 divided by wave length, is introduced in place of w 


7 


because gust disturbances are essentially space disturb- 
ances rather than disturbances in time. The several 
measurements that have been made to establish this 
spectrum will be reviewed later. The second sketch 
7°(Q2) represents the amplitude square of a specified 
airplane response, such as the airplane normal accelera- 
tion, to sinusoidal gusts of unit amplitude and of fre- 
quency 2. The sketch shown is a typical normal ac- 
celeration response for a moderately flexible airplane. 
Airplane modes of response in this representation show 
up as peaks such as the free body and fundamental wing 
bending modes illustrated. The output spectrum %)(Q) 
is obtained in accordance with Eq. (S) as the product of 
the first two curves ®,(Q2)7°(Q) and gives, for example, 
the spectrum of normal acceleration or a spectrum of 
stress. The significance of this spectrum in regard to 
the characteristics of the airplane response is considered 


in the next section. 


Relations Between Spectrum and Time History 


In addition to the spectrum, specific statistical char- 
acteristics of the response time history are usually of in 
terest. Thus, for example, the largest value or the num 
ber of exceedances per second of given intensity levels 
are frequently required. In some cases, the spectrum 
also provides information on these statistical charac- 
teristics and it thus appears appropriate to review the 
information that the spectrum may provide in this re- 
spect. 

It will be recalled, Eq. (7), that the integral of the 
spectrum determines the mean square value. The root 
mean square value 0 = V y2(7) may thus be obtained 
from the spectrum and provides a useful linear measure 
of the disturbance intensity. If the disturbance y(/) 
has a normal or Gaussian probability distribution with 


zero mean, the probability density p(y) is given by 

p(y) = (1/ovV/2r)e" se (9) 
The probability distribution defines the proportion of 
total time that the disturbance is at any value and is in 
this case completely specified by the root mean square 
value o. Fung, in reference 9, has also pointed out that 
for the Gaussian case, the largest value of y for a given 
time is directly proportional to the root mean square 
value. 

Additional relations between the power spectrum and 
the time history characteristics have been investigated 
in an extensive study by S. O. Rice.6 For the special 
case in which the disturbance is completely Guassian 
(one in which the joint probability of y(t), y(t + 7), 
y(t + re), ete., is Gaussian), Rice has derived a number 
of relations which appear useful in aeronautical applica- 
tions and are particularly significant for fatigue studies. 
Of these, the more important appear to be expressions 
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Fic. 1. Gust response determination 


for the average number of times per second, V,(y), that 
a given value of y is crossed with positive slope and the 
average number of peak values (maxima) per second, 
N,(y), that are above a given value of y. The average 
number of times per second that a given positive value of 
y is crossed with positive slope (or a given negative 
value is crossed with negative slope) is given by 


) 


NAy) = (1/27)(01/o)e or (10 


where o; is the root mean square value of y’(f) and is re- 
lated to the spectrum of y(‘) by the relation 


on = Vy’)? = | w"P(w) dw (10a) 


The expression for the number of peaks per second that 
are above a given value y (or the average number of 
minima or ‘‘valleys’’ per second that are below —y) is 


somewhat complicated but for large values of y 
N,(y) = NAY) (11) 


For most practical cases, it appears that this approxima- 
tion is quite good for y > 2c. The total number of 
peak (maxima) values of y (or the number of minima) 
which includes the maxima for positive and for negative 


values of y, is given simply by 
N, = (1/22)(o2/a;) (12) 


where 


02 = 7 [y"()]? = | wid(w) dw (12a) 
Jo 
It may be remarked that, at least for some cases, the 
loads experienced by an airplane appear to follow a 
simple Gaussian distribution. The foregoing relations 
between the spectrum and the peaks, Eqs. (10) to (12), 
which apply to completely Gaussian disturbances, have 
been tested empirically in several cases with some suc- 
cess and will be used in the present applications. 
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Estimation of Spectra from Experimental Data 

General Considerations..-As a final aspect of the 
theory, consideration should also be given to the prob- 
lem of spectra estimation from experimental data. The 
expressions for the power spectrum, Eqs. (3) or (6), in 
volved infinite integrals. Experimental determinations, 
on the other hand, must of necessity be based on finite 
and, in practice, reasonably short record lengths, 
thereby introducing a degree of approximation. In ad 
dition, evaluations are generally restricted to uni 
formly spaced values of time. The problems introduce | 
by these finite approximations are not trivial and some 
of them have apparently not been appreciated until 
recent years. 

Spectra determinations can, in practice, be made 
through numerical procedures or by the use of analog 
equipment. The numerical methods require extensive 
calculations; the analog methods are quicker, but not 
as precise. The following discussion is given in terms 
of numerical evaluations although similar considerations 
are involved in the use of analog methods. It is based 
on the results obtained by John Tukey and his col- 
leagues at the Bell Tellephone Laboratories, who by 
1949 had developed a mathematical foundation for the 
technique of deriving and interpreting spectra esti 
mates.* Because these results are not in readily ac 
cessible sources, this summary of their significance and 
method of application is included in the hope that it 
may prove useful in aeronautical applications. 

Two fundamental problems are introduced by the 
aforementioned finite approximations and concern the 
frequency resolution and the statistical reliability of the 
estimates. 

Frequency resolution: Frequency resolution 1s con 
cerned with the ability to accurately distinguish between 
power at a given frequency and power at adjacent fre- 
quencies. It may be shown that, because of the finite 
record length and the numerical procedures involved, 
some of the power at a given frequency has a tendency 
to ‘‘diffuse’’ or show up as power at other frequencies. 
As a consequence, estimates of power are at best 
weighted averages of the power over a finite frequency 
range. Also, if the power is to be estimated over the 
frequency range 0 to say w , power present at higher 
frequencies than wo will tend to show up as power at 
frequencies below wo. Specifically, power in the vicinity 
of 2a) + w, 4a + w, 609 + w, etc., Shows up as power at 
w. Tukey has denoted this characteristic by the term 
“aliasing.”’ 

Statistical reliability: This problem arises primarily 
from the finite length of the record since two adjacent 
sections of a record, like two series of coin tosses, will 
generally not yield identical results. The statistical 
reliability of the estimates must be established in order 
to discriminate between real and statistical fluctuations. 

The effects of the finite approximations on measured 
spectra may be interpreted as equivalent to convoluting 
the true spectrum with a unit area “‘band-pass filter’’ as 


given by the following relation 


rs 8) (w) = 
I meas, \“/ 


P,, te (W@W A (Ww Ww) dw 


where A(w) is the characteristic weighting function ot 
the so-called filter. Thus, the measured spectrum is a 
“smeared”’ representation of the true spectrum, though 
the total power remains unaltered. The actual char 

icteristics of the filter A(w) depend upon the precis¢ 
numerical procedures used and can to some extent bx 
modified. It is desirable, therefore, to tailor th 

numerical procedures so as to yield an effective filter 
which minimizes both the tendencies for power di 

fusion and the statistical fluctuations of the estimate 

On the basis of these considerations Tukey? arrived at 
a simple and quite practical numerical calculation 
procedure for spectra estimation. The application « 

the calculations procedure for a given case requires 
the choice of the following quantities: 

(a) T record length in time 


(b) At = T/n 
(c) m = number of uniformly spaced points 


time interval of record reading 


over the frequency range of in 
terest at which power estimates 


are desired 


Criteria for the appropriate choice of values for these 
quantities in applications will be indicated in the subse 
quent discussion which in the interest of concreteness is 
deferred until after the calculation procedure is outlined 
and the characteristics of the estimates are described. 
The starting point for the development of the Tukey 
procedure is the communication sampling theorem 
which says that sampling at intervals of time Af vields 
all the information obtainable on frequencies up to w 

m/ At or frequencies which are sampled at least twice per 
wave length. The actual calculation involves the steps 
outlined below. 

Tukey's Calculation Procedure for Spectra Estima 
tion. -(1) Determine the auto or serial correlation co 
efficients in accordance with Eq. (5b) from the succes 
sive values y\, Vo... ¥, spaced at uniform intervals A, 
by means of the following expression 


R, = . 2 VeVe+ 3 (f= © 1,2...) 


where 
R. = Rlr = rf) 


It will be noted that only the truncated autocorrelation 
function covering 7 = 0 to 7 = mAt is used. 

(2) Initial or “‘raw”’ estimates of the power density 
are first determined by the Fourier series cosine expan- 
sion in accordance with Eq. (6) 


At m— } ” 
L, = Ro +2 > R, cos _ bed + R,, cosrr (14 
7 q l m 


where 
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L, = Lie = (re/nAt) 


These estimates, however, have an effective filter which 
has the undesirable character of appreciable side band 
areas; thus, a wide diffusion of power results. 

(3) Final or “‘smoothed”’ estimates of power density 
which are estimates based on a more desirable and 
sharper effective filter are then determined by 


&) = 0.5415 + 0.46L, 
®, = 0.23L,_, + 0.541, + 0.23L,4, $ (15) 
$, = 0.46L,,—, + 0.542 


Frequency Resolution. The values of ®, obtained in 
this manner consist of m + 1 power estimates uni- 
formly spaced over the frequency range 0 to m/At. 
These estimates can be considered as estimates of the 
weighted average power over the frequency bands 
(r — 1)ar/mAt and (r + 1)r/mAt. The actual effective 
filter shape or weighting function is roughly triangular 
with a base equal to four times the time between power 
estimates (42/mAt). Because of the large overlapping 
of the weighting functions for adjacent estimates, adja- 
cent estimates are not statistically independent; 
however, every other power estimate, for example, ®», 
,, etc., can be considered statistically independent. 
Thus the effective frequency resolution is twice the 
band width between estimates or 27/mAt. It will be 
noted that increasing the number of separate power 
estimates, m, narrows the effective filter band width 
and thus has the effective of increasing the frequency 
resolution. The price of this increased resolution as will 
be seen, is, however, a loss in statistical reliability for 
the individual estimates. 

Statistical Reliability. For the case of a Gaussian 
disturbance, Tukey has shown that estimates of ® de- 
rived in the foregoing manner scatter about the true 
value and have a frequency distribution that is ap- 
proximated by the so-called ‘“‘chi-square”’ distribution." 
For spectra that are known to be relatively flat, the re- 
liability of measured values of power is dependent only 
upon the sample length and the number of separate esti- 
mates of power made. The reliability of each estimate 
may be expressed in terms of the quantity, k, the so- 
called statistical number of degrees of freedom, where 


k is given by 
k = (n — n/5)/(m/2) (16) 


The reliability of estimates may be indicated by con- 
fidence bands. These bands provide a measure of the 
range within which the true value of the power lies for 
a given probability level. Fig. 2 gives the 90 and 95 
per cent confidence limits of the quantity ®,/®,, , where 
?,. represents the average of the ‘true’ power over the 
band width, as a function of the number of degrees of 
freedom, k. Thus, for example, for k = 100, the follow- 
ing 95 per cent confidence limits for &,, are obtained: 
0.77Pmeas < Pay. < 1.398 


meas meas 


That is, the true average power can be expected to be 
between 0.77®Pyyeas and 1.59,, with a probability of 
95 per cent. To obtain greater reliability, more degrees 
of freedom are of course required. This may be ob 
tained by either larger record length or at the sacrifice of 
frequency resolution, by averaging over wider frequency 
bands (smaller m). 

For spectra that are not flat, the number of degrees of 
freedom for some of the estimates may be reduced from 
that given by Eq. (16). The reduction in each case 
depends upon the shape of the spectrum over the band 
width for which the estimate applies. A useful rule of 
thumb given in reference S for determining the effective 
number of degrees of freedom for each of the estimates, 
for this case is 

(A; + Ao +... Ax)? ” 
si 3 : (17) 

Ae tt Ae Tote 
where the A's represent the true power at the 2”, m in- 
termediate frequencies uniformly spaced over the fre- 
quency band width for which the estimate applied, 
(yr — 1)w/mAt to (r + 1)r/mAt. If the X's are equal in a 
given frequency band with, Eq. (17) yields 2”/m for the 
number of degrees of freedom which is close to the result 
given by Eq. (16). If all the \’s except one are small (a 
narrow peaked spectrum), then Eq. (17) yields a value 
of 1, thus leading to an extremely poor power estimate. 
The reduced number of degrees of freedom obtained 
from Eq. (17) may be used with Fig. 2 to obtain confi- 
dence bands as already indicated. In practice, the true 
spectrum is frequently not known; therefore, the 
measured power estimates have to be used as a guide. 
It is of interest to note that the reliability of the 
measured values of total power, as contrasted to the re- 
liability of the individual powers, can also be established 
with Eq. (17) by letting the \’s represent the true power 
at » + | uniformly spaced frequencies from 0 to the 
highest resolved frequency 7/A/, or more practically by 
using the (m/2) + 1 independent power estimates 

P,, bs, etc., as a guide to the values of X. 

Choice of At, m, and n._-The foregoing results in ad- 
dition to providing estimates of reliability, also provide 
a systematic basis for designing experiments for spectra 
determination. Through their application, experi- 
ments and calculations may be tailored to achieve both 
desired reliability and frequency resolution by the ap- 
propriate choice of the values of Af, m, and n. The 
following criteria are indicated as a guide to spectra 
estimation. It will be recalled that the highest resolved 
frequency is given by mw/Aft. Thus, the highest fre- 
quency desired determines the reading interval which 1s 
given by 

At = /wo (1S) 
In practice, this value should be chosen to assure that 
the power at the higher frequency is negligible since 
otherwise this power would tend to leak into the power 


estimates at the lower frequency. 
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STATISTICAL DEGREES OF FREEDOM,k 


Fic. 2. Confidence limits for &,/®,,. as a function of k (based on 
chi-square distribution ) 


The appropriate choices of m and » depend on both 
the desired frequency resolution and the statistical re- 
liability, and are intimately related to the spectrum 
shape. If the spectrum is known to be relatively flat, a 
small number of power estimates m will be sufficient to 
describe the spectrum. The appropriate value of ” can 
then be obtained from Eq. (16) or is roughly km/2 and 
depends on the statistical reliability desired as given by 
Fig. 2. For spectrums that are not flat, the value of 1 
should be increased in accordance with the reduction of 
the number of degrees of freedom as indicated by Eq. 
(17). If the spectrum has a sharp peak and a reliable 
measure of the peak power is required, the width of the 
peak will determine the appropriate frequency band 
width for the power estimates and thus the appropriate 
value for m; that is, the interval 27/mAft should 
correspond roughly to the spectrum width near the top 
of the peak. The appropriate value of m is then ob- 
tained from Eq. (16) or (17), where the value for k is 
chosen from Fig. 2 for the level of reliability desired. 


AVAILABLE MEASUREMENTS OF ATMOSPHERIC 
TURBULENCE 


In order to apply the foregoing techniques to the cal- 
culation of airplane response to continuous rough air, it 
is necessary to define the spectral characteristics of at- 
mospheric turbulence. For the purposes of the present 
applications, the available data on the spectrum of at- 
mospheric turbulence will be briefly reviewed in order to 
provide a guide to the selection of a representative tur- 
bulence spectrum. The information on the spectrum of 
atmospheric turbulence is quite limited. The first 
measurements of the turbulence spectrum in the free 
atmosphere were made by Clementson in 1950.' A 
few subsequent studies have been made and include 
some results obtained by the Douglas Aircraft Co."! 
and unpublished NACA results. The results obtained 
in these several studies are summarized in Fig. 3. A 
log-log plot is used in order to cover the broad range of 
frequency and power density represented. 


The spectra shown represent absolute values, no 
attempt having been made to normalize the reported 
results to the same mean intensity level. In some cases 
the results shown are for the average of several test runs 
and cover several different weather conditions. All but 
one of the spectra are for vertical gust velocity as meas 
ured with an airplane in horizontal flight. The one ex 
ception, denoted by /7/, was obtained from the high fre- 
quency fluctuations of a sensitive air-speed system and 
thus represents the spectrum of horizontal (or longi 
tudinal) gust velocity. Although the mean intensity 
level differs considerably from one test to another, as 
judged by the height of the curves, there is reason to 
believe that the conditions represented by these results 
do not cover the more severe conditions of turbulence 
such as found in thunderstorms. In fact, for a severe 
thunderstorm, it might be expected that the turbulence 
would have a power level perhaps ten times as severe 
as those shown. 

Perhaps the most significant point to note about the 
curves is that in most cases the power decreases at a 
uniform rate with increasing frequency. For the range of 
frequencies represented 0.002 < 2 < 0.5, the indications 
are that the spectrum may be reasonably approximated 
by a c/Q? type curve, c being a constant which evidently 
depends on the intensity level. It is of interest to note 
that available results derived from the theory of iso 
tropic turbulence and wind tunnel measurements of 
turbulence yield spectrum shapes at the higher fre- 
quencies that conform with the 1/Q2° shape indicated 
herein. In view of these indications, a spectrum shape 
of 1/Q° will be assumed for the range of frequencies 
0.002 < 2 < 0.5 (gust wave lengths of roughly 10 to 
3,000 ft.) and this spectrum will be used in the subse- 





o NACA 

4 NACA 

© NACA (h) 
6 GLEMENTSO 
---DOUGLAS 











rs 





% 
1000 ‘ 
\ 








100 
(FT/SEC)* x 
-RAD/FT ’ a 


lo + 





® 









































| | 
OO! Ol A! 10 
2, RAD./F 


Fic. 3. Measured power spectra of atmospheric turbulence 
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Fic. 4. Effect of c.g. position 


quent calculations. For each of the calculations, ap- 
propriate values for the turbulence intensity will be in- 
cluded as available. 


APPLICATIONS TO Gust LOADS ON AIRPLANES 


In this section application of the foregoing theory is 
made to investigate three specific gust problems. The 
cases to be treated are concerned with the effects on gust 
loads and stresses of: (1) center of gravity position, (2) 
low damping in short period, (3) wing bending flexibil- 
ity. 


Effect of Center of Gravity Position on Gust Loads 


This application deals with the effects of airplane 
center-of-gravity positions on the loads encountered in 
continuous rough air. Several years ago, the NACA 
conducted some side-by-side flight tests using two jet 
fighter aircraft differing only in center-of-gravity posi- 
tion. The overall test results for an air speed of 450 
m.p.h. reported in reference 12 indicated roughly a | 
per cent increase in center-of-gravity normal acceleration 
per | per cent rearward movement of the airplane center 
of gravity. In order to determine whether these experi- 
mental results agreed with what would be expected 
from theoretical considerations, calculations were made 
for the two test conditions using the power spectrum 
approach. 

The airplane equations of motion for the stick fixed 
condition were solved and the normal acceleration fre- 
quency response function for gusts determined. For 
these calculations, the airplane was assumed rigid but 
free to move vertically and pitch. The amplitude 
squared of the frequency response functions for normal 
acceleration are shown in Fig.4a. The curves shown are 
for the actual position of the accelerometers used in the 
tests which because of space limitations were located 
several feet forward of the airplane center of gravity. 
The peaks in both cases correspond essentially to the 
airplane short period mode. It will be noted that as the 
center of gravity moves rearward, the peak frequency 


response moves to a lower frequency and decreases 
somewhat. 

Using the 1/Q2* gust spectrum previously indicated 
and a rough estimate of the average turbulence in- 
tensity, the power spectrum of airplane normal ac- 
celeration for the two center of gravity positions were 
obtained and are shown in Fig. 4b. As would be ex- 
pected, the peak of the response spectrum moves to 
lower frequency as the center of gravity position moved 
rearward. However, due to the shape of the gust spec- 
trum, the acceleration spectrum for the rearward center 
of gravity position is now larger and covers more area. 
In fact, the ratio of the root mean square values for the 
rearward and forward center of gravity positions is 1.13 
indicating a 13 per cent increase as the center of gravity 
moved rearward. This increase is in substantial overall 
agreement with the test results of reference 12 which in- 
dicated a roughly 10 per cent increase for these con 
ditions. 

In order to obtain a more detailed comparison be- 
tween calculated and test results, the previously men- 
tioned relations between spectrum and number of peaks 
derived by Rice [Eqs. (10) to (12)| were applied to the 
spectrums of Fig. 4b and the average number of peak 
accelerations per second that exceed given values were 
determined. These calculated results are shown in Fig. 
5) and indicate small difference at low values of An in- 
creasing to around 15 per cent at lg. The comparable 
results from the flight tests data also shown in the figure 
indicate essentially similar trends. The calculated re- 
sults in this case thus appear to closely reflect the mag- 
nitude and character of the load’s increase as the center 
of gravity position is moved rearward. 


Effects of Short Period Damping 


The second application to be described concerns the 
effects on gust loads of airplane dynamic stability, par- 
ticularly the effects of low damping in short period. 
Recent theoretical studies,” ‘ have indicated that low 
short period damping might seriously amplify the gust 
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loads in continuous rough air. During the last year 
some test data bearing on this problem have been ob- 
tained by the NACA. A 45-degree swept-wing tailless 
missile was flown through rough air at low altitudes and 
high speeds and the results indicated that a very sub- 
stantial amplification in gust loads was present. 

The dynamic stability characteristics of the missile 
are summarized in Fig. 6. The period, shown as a func- 
tion of air speed, is seen to be very short (0.11 sec.) and 
decreases somewhat with speed. The time to damp to 
one-half amplitude is about 0.15 sec. at the lower speeds 
covered but increased abruptly in a small speed range to 
a value of about 0.24. The vertical dashed lines indi- 
cate three of the test air-speed conditions |, Io, and V3, 
which covered a 10 per cent variation in air speed in the 
region of damping deterioration. The results obtained 
for these speeds will be described in some detail. 

For comparison with the test results, calculations of 
the expected spectrums of normal acceleration were 
made for these three air speeds and are shown in Fig. 
7a. These were obtained from calculated estimates of 
the missile normal acceleration frequency response func- 
tion and the gust spectrum previously indicated. In 
these cases, the intensity of the turbulence was esti- 
mated from airplane survey flights over the missile 
firing range just prior to the tests. These spectrums 
show a progressive movement to higher frequency as the 
air speed increases, directly reflecting the variations in 
the short period frequency. The spectrums are also 
progressively more peaked, reflecting the deterioration 
in the damping of the short period frequency or the re- 
sultant narrowing band pass character of the missile 
response. The relative load levels for the three air- 
speed conditions as measured by the root mean square 
values are also given and show a relatively large in- 
crease, roughly from 0.4g to 0.65g, for the 10 per cent 
increase in air speed. 

For comparison, the measured power spectrums of 
normal accelerations for the same three test air speeds 
are shown in Fig. 7c. These spectrums were obtained 
by applying the Tukey technique, Eqs. (13) to (15), to 
the measured time histories. In each case, only very 
limited samples were available (roughly 2 sec.), with 
the result that only limited statistical reliability could 
be obtained. Nevertheless, the calculated and measured 
root mean square values are seen to be in reasonable 
agreement. In obtaining the power spectrum esti- 
mates, statistical reliability considerations dictated a 
choice of a rather wide frequency band width; each 
point represents a weighted average over roughly 3.2 
cycles per second (+1.6 cycles per second about the 
values shown). The estimates may therefore be ex- 
pected to provide a “‘smeared”’ picture of the real spec- 
trum, and this is borne out by the fact that measured 
results are flatter in character than are the calculated 
results. 

In order to make a more fair comparison between the 
calculated and measured results, the calculated results 


AERONAUTICAL 


SCIENCES JANUARY, 1953 

of Fig. 7a were modified by the application of a filter 
equivalent to that present in the measured estimates. A 
triangular filter A(/) of 3.2 cycles per second base width 
and unit area was used and the modified results ob 
tained by the relation: 


. 
P04. (f) = / Poy (SIAC — fi) dhi 
Jo 
The results obtained with this modification are shown 
in Fig. 7b and are seen to be in good overall agreement 
with the measured spectrums of Fig. 7c. In fact, con 
sidering the limited statistical reliability present, the 
agreement seems much better than might have been ex 
pected. As mentioned previously the root mean square 
values are not changed by this modification. 

The overall results of the missile test and calculations 
are summarized in Fig. S which shows the variation of 
This 


solid curve represents the calculated results, the circles 


root mean square acceleration with air speed. 


represent the measured values obtained for eight test 
air speeds. The experimental results show considerable 
scatter as would be expected from the limited test data. 
Nevertheless the experimental data follow the general 
trends of the calculated curve. The increase in root 
mean square value for the air-speed range covered is 
roughly 85 per cent, far in excess of the 30 per cent in- 
crease that would be expected from the direct effects of 
the increase in air speed and the changes in the slope of 
lift curve. The additional 55 per cent increase in root 
mean square value appears to be a direct consequence 
of the deterioration of the short period damping in the 
test range. These results thus appear to confirm the 
theoretical indications of large load amplification asso- 


ciated with low damping. 


Effects of Wing Flexibility 

The third and final application considered herein is 
concerned with the magnification in gust load stresses 
that results from wing flexibility. Flight tests had been 
made in clear rough air with the three aircraft pictured 
in Fig. 9 to establish what the effect of wing flexibility is 
in practical cases, and it was of interest to see how well 
“power spectrum” calculations could predict the flight 
test results.’ Relatively, the three aircraft employed 
were judged to be respectively rather stiff, moderately 
flexible, and somewhat flexible. The rough air flight 
tests confirmed this comparison and indicated that bend- 
ing flexibility was the chief wing flexibility of concern. 
Specifically, the flight tests indicated that for the three 
aircraft the incremental root bending stresses due to 
gusts were, respectively, 5, 20, and 28 per cent greater 
on the average than what the stresses would be if the 
airplanes were rigid. 

In this application calculated results were obtained 
by using the atmospheric spectrum given in reference 2 
as the input. The frequency response functions for 
bending stress were obtained by reducing the analysis 
of reference 13, which involves the two degrees of free 
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dom, vertical motion and wing bending, to the special 
case of a continuous sinusoidal gust disturbance. The 
response functions were evaluated for weight conditions 
and spanwise stations representative of those used in the 
flight studies and are shown in Fig. 10. The solid curve 
is for the flexible aircraft; the dashed curve, for the air- 
craft considered rigid. The curves show quite clearly 
the different bias that each airplane has toward various 
frequency components of the atmospheric motions. 
The first hump is associated with vertical translation of 
the aircraft; the second hump, with wing bending. 
The spectra for bending stress response, obtained by 
multiplying the frequency response curves by the input 
spectrum, of course show that the curve for the flexible 
case overshoots the curve for the rigid case by an 
amount consistent with the frequency response curves. 
The area of this overshoot is a direct measure of the am- 
plification in mean square bending stress that results 
from wing flexibility. 

To obtain an amplification or flexibility measure 
more directly comparable to the values obtained from 
the flight tests results, the following procedure was 
used. The number of peaks per second that exceed a 
given stress level were calculated in accordance with 
Eqs. (10) and (11) for both the flexible and rigid con- 
ditions. At given numbers of peaks per second, the 
ratio of the stress for the flexible airplane to the stress 
for the rigid case was then formed. This ratio was used 
as a measure of the amplification due to flexibility. The 
ratio varies somewhat with stress level, but for com- 
parison with flight tests results the ratio was chosen at 
a stress level equal to about twice the root mean square 
stress. 

Fig. 11 shows the correlations of the results obtained 
by the harmonic-analysis approach with flight results. 
The ordinate is the amplification values for stress that 
were obtained in flight, the abscissa is the ratio ob- 
tained from the harmonic-analysis theory as explained 
in the previous paragraph. The good correlation 
shown by this plot is, as with the previous applications, 


very gratifying. 


CONCLUDING REMARKS 


The foregoing discussion has served to indicate that a 
fairly complete mathematical technique now exists for 
the treatment of random type disturbances such as 
occur, for example, in the atmosphere and perhaps in 
buffering fluctuations. The technique includes: 

(1) The mathematical representation of the random 
disturbance by means of the power spectral density 
function. 

(2) The use of this function in conjunction with the 
response characteristics of a system for determining the 
power spectrum of the system response to this dis- 
turbance. 

(3) The relation of the power spectrum to specific 
statistical characteristics of the response such as the 
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mean square value, and in particular cases, the proba- 
bility distribution and the distribution of peak (maxi- 
mum) values. 

(4) Methods both for 
measurements and for establishing the reliability of such 


estimating spectra from 
estimates. 

As a preliminary to the application of the foregoing 
techniques to some gust load problems, a brief review of 
the available information on the spectrum of atmos 
pheric turbulence was presented. The available data 
indicated considerable variation in the intensity of the 
turbulence spectra reflecting a variation in turbulence 
intensity among the various weather conditions. The 
spectra in all cases, however, appear to indicate a uni 
form rate of decrease in spectral density as the frequency 
increases (gust wave length decreases). For the range 
of wave lengths covered, 10 to 3,000-ft. wave lengths, 
the spectra appear to approximate a shape of 1/2’, 
where 2 is a reduced frequency in radians per foot. 
This spectrum shape appears in conformity with avail- 
able results in the theory of isotropic turbulence and 
wind tunnel measurements of turbulence spectra. 
This spectrum shape is used in subsequent gust load 
calculations. 

Several applications are then presented in order to 
indicate some effects on gust loads and stresses of: (1) 
airplane center of gravity position, (2) low damping in 
short period, and (3) wing bending flexibility. In each 
case, experimental confirmation is also given. 
cases considered, the effect of a rearward movement of 
the center of gravity is shown to result in about a 15 


For the 


per cent increase in the larger loads in rough air. The 
deterioration of short period damping is shown in a 
particular case to result in a 55 per cent amplification of 
the loads. The effects of wing bending flexibility in 
stresses are shown to be small, 5 per cent, in the case of 
one particular airplane, but to yield sizable stress ampli- 
fication, 20 and 28 per cent, respectively, for two other 
airplanes. From the rather good agreement that was 
obtained between the experimental and theoretical re- 
sults, it appears that the generalized harmonic analysis 
approach can be an effective means for extending the 


art of gust load analysis. 
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ABSTRACT 


Phe reflection and diffraction of strong shocks around corners of 
arbitrary, finite angle, and the resuiting pressure and density 
fields have been calculated by two different methods 

The second method, which is a hyperbolic procedure, appears to 
be useful in various types of nonlinear problems. At the end a 
modified procedure of the second method is outlined which prom- 


ises faster convergence 


INTRODUCTION 


HE REFLECTION AND DIFFRACTION of strong shocks 
yee concave corners of arbitrary, finite angle, 
and the resulting pressure and density fields have been 
calculated by two different methods. 

In the first procedure, the problem is formulated in 
conical coordinates; as a consequence, the nonlinear 
equations of motion become essentially elliptic, and 
boundary conditions have to be imposed all along the 
boundary of the domain of flow disturbance. Some 
conditions are “‘floating’’ (on shocks), others are on 
fixed boundaries, and their number is different on dif- 
ferent boundary sections. The solution has been ob- 
tained in numerical form, as well as by a series expan- 
sion about a singular point, which appears to be the 
key to the understanding of the pattern. Agreement 
with shock-tube interferograms is satisfactory. 

In the second procedure, the problem is formulated 
in physical coordinates x, y, ¢; consequently, the equa- 
tions of motion are hyperbolic, and only one boundary 
Now, there are 
no ‘floating’ conditions whose locations would involve 
trial-and-error Shocks and _ slipstream 
come out as part of the general solution, in form of steep 


procedures. 


gradients of density (and pressure) at the correct loca- 
tion. At present, this problem is being coded for the 
IBM ‘Defense Computer.”’ 

The second method appears to be useful in several 


other types of nonlinear problems. 


ANALYTIC FORMULATION OF THE PROBLEM: 
ELuiptTic METHOD 


Fig. | represents the diffraction of a strong shock (of 
pressure ratio po/p; = 0.36) around a concave corner! * 


Presented at the Aerodynamics Session, Twenty-Second Annual 
Meeting, IAS, New York, January 25-29, 1954 

* Investigation sponsored by the Mechanics Branch of the 
Office of Naval Research. 

+ Associate Professor of Aeronautical Engineering 

t Research Associate in Aeronautical Engineering 
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Both 
experiment and theory show that for strong shocks and 


of @ = 38° and the resulting density field p/p». 


finite corners, as indicated, Mach reflection should re- 
sult. 
in the problem, the p- and p-field will grow with time, 


Since no characteristic, physical length appears 


geometrically similar to itself; hence it is appropriate to 
introduce conical coordinates £ = x/t, 7 = y/t, in which 
While 


glancing incidence could be solved by linearized theo 


the pattern appears static. conditions at 
ries,* the Mach reflection of shocks from finite corners 
leads to a difficult nonlinear boundary value problem, 


which may be formulated as follows: 


a) Differential Equations‘ 


Conservation of mass: 


(pU); + (pV), + 29 = 0 (1 
Conservation of momentum: 
U-U,;+ V-U, + U = —p,/p (2 
U-V;+ V-V,+ V= p,/p (3 
Conservation of energy: 
U-S; + V-S, = 0 (4 


¢ 
< 


Here p, p, and S are pressure, density, and entropy; 
furthermore U) = u — and V = v — 7 are the com- 
ponents of a vector VY and wu, and v are the actual fluid 
velocities at any point &, 7. 

From Eq. (4) it follows that the ‘‘stream lines” 

V/U = (dn/dé)s (5 
are isentrops. In accordance with the definition of 
U and I’, the isentrops are straight lines pointing radi- 
ally, toward the corner 0 in domain (0) in front, and 
pointing toward C in domain (1) behind the incident 
shock (see Fig. 1); here C is at a distance u,/a, from 0, 
where 1,/a,; is the Mach Number of the undisturbed 
flow in domain (1). After substituting the entropy S in 
terms of p and p, the above system represents four 
equations for the four unknown dependent variables 
Ul’, V, p, and p. Noting that Eqs. (1) to (4) differ in 
form from the stationary equations only by nondiffer- 
ential terms, it is seen that two of the four character- 
istics are a real double root—-namely, the isentrops 
the other two are imaginary. 
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(b) Boundary Conditions 


The location of the reflected and Mach shocks and 
that of the slipstream are not known a priori. At these 
boundaries of unknown location, the four conservation 
laws must be fulfilled. This yields four conditions 
on the shocks and two conditions on the slipstream, the 
other two being automatically fulfilled. On the sym- 
metry line and on the wall only one condition is to be 
namely, the usual requirement of tangential 
flow velocity. To solve the problem, a shock and slip- 
stream configuration has to be assumed. Actually, 
the correct location of these discontinuity surfaces 
The problem would 


imposed 


follows from the desired solution. 
be overdetermined if the fulfillment of all four conserva- 
tion laws would be required at the selected boundary. 
Consequently, one boundary condition has to be 
dropped temporarily; it will serve later on as a check 
for the correctness of the assumed boundary form. In 
other words, a_ trial-and-error procedure becomes 
necessary, which, however, can be reduced by taking a 
hint from the experimental configuration. Taking the 
preceding considerations into account, there must be 
three conditions on the shocks, one on the wall, and one 


on the slipstream. 
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To carry out a numerical solution of this nonlinear 
system, in difference form, the boundary values of all 
four dependent variables have to be known along the 
whole boundary. (Otherwise, there would not be 
enough data for the determination of the interior net 
points.) Hence four boundary conditions are needed 
everywhere, and the aforementioned true boundary con- 
ditions have to be supplemented by so-called dependent 
conditions which are an expression of the validity of the 
differential equations at the boundary. 

Hence, on the reflected and Mach shocks, the three 


true boundary conditions may be written: 


(On) B = Ag, (Ens Nshy Po, ly Poy 1) | 
(Q,) B => Bo, (Esa Nsn) ; (6) 
(p)pB sa Co, (Esny Ash» Po ly Puy 1) 

where subscripts B = boundary, sh = shock, 0 = in 


front of and 1 = behind incident shock; s is parallel and 
n is perpendicular to the boundary. Furthermore, the 
functions A, B, C and the quantities on the right are 
known. 

A linear combination of the momentum equations 
yields the following one dependent boundary condition 


for the boundary value pz: 





Solid lines: Isentrops 
Dashed lines: Isopycnics 
Po/P, = 36 
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oP) 2 Bi . oU 
= - COs~ ma l ] -lg sh 
( 4 p COS® Qs, \| a ga (5 
OV . 4, OU oV')) . 
+ [U-tg ay, — | —_— > 
s a IE On i} 


where tg ay, = (dn/dé),,. Since, in difference form 


(Op/On)p =) (pe — pr)/An (7a) 


px may be expressed by Eq. (7) in terms of the pressure 
in the interior p, and other quantities appearing any- 
how in the equation at J. 
On wall and slipstream, the one true boundary condi- 
tion is 
V/U = tang 6(&, n) = (dn/dé)s (S) 


/ 
where @ is the angle of the boundary with the é-axis. 

In accordance with Eq. (8), for both boundaries Q,, = 
0, by definition. This proves that wall and slipstream 
are isentrops. 

The three dependent conditions at the wall are 


Op/on = 0 

Op — pu 4 oF) sec" 6 (9) 
Os Os Os 

p/p" =o ) 


On the slipstream, the three analogous conditions are 


Op/On = pl” sec® 6-(db/ds) 


op = —pl’- seca + Ve+e| (9a) 


Os Os Os 

p/p’ =r 
here tg 6 = (dn/dé);_; indicates the direction of the 
slipstream. Furthermore, o; and o2 are the two respec- 


tive entropy constants. 


GENERAL DISCUSSION OF THE SOLUTION 


The preceding formulation of the problem is sufficient 
to carry through a numerical difference scheme (to be 
described elsewhere) if a plausible boundary configu- 
ration is assumed. It is evident that the system of the 
four unknowns U, V, p, and p is uniquely determined. 
At each interior net point, Eqs. (1) to (4) have to be 
satisfied; at each boundary point, the four true and 
dependent conditions are to be fulfilled. The solution 
was obtained in two steps. First, the Q-lines or isen- 
trops were computed by using the experimental iso- 
pyenics. Eqs. (1) to (4) were reduced to two equa- 
tions for U and V, by eliminating p from Eqs. (2) and 
(3) and considering the p-field as given. This yields a 
valuable first guess. Afterwards, the p-field was con- 
sidered as unknown, and Eqs. (1) to (4) were used to 
determine all four unknowns. No discrepancy between 
the computed and experimental isopycnics could be 
discerned, within the accuracy of the calculation; nor 


was there any deviation between theoretical isobars 
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ANALYTIC RESULTS NEAR S AND O 


Fics. 2 and 3. 


and isopyenics, except in the vicinity of S, which is 
discussed in Figs. 2 and 3. The key to the under- 
standing of the disturbance field NOST is the behavior 
of the isentrops near S, which shall be analyzed in the 
following: 

(A) By Eq. (8), the slipstream is an isentrop; more 
exactly two isentrops because of different values for the 
entropy constant o on either side. The symmetry line 
and wall section NOS represent a third, and wall sec- 
tion MS a fourth isentrop; all of which meet at the 
point S on the wall. The value of ¢ on all four isen- 
trops can be seen to differ, due to different strengths 
of the shocks at NV, 7, 7’, and M. Outside the dis- 
turbed region NOST, the streamlines Q are rays con- 
verging toward C. In the linearized case of a shallow 
corner ¢, in which the disturbance is of O(€), these rays 
continue radially all the way through. For large cor- 
ners the radial character of the isentrops will be dis- 
turbed. Since the resulting curves cannot turn around 
and continue outwardly, they have to meet all at the 
intersection point of slipstream and wall. In other 
words, at S, there is a directional singularity of the 
entropy. 

(B) From the ‘‘ray behavior’ of the isentrops in the 
vicinity of S, it follows that Q = 0 for any isentrop at 
S; otherwise, there would be a Q,-component perpendi- 
cular to the wall. Hence the Q-values along any isen- 
trop increase from 0 to some value smaller than 1, on the 
reflected shock. As will be shown later, it is a necessary 
consequence of Eqs. (2) and (3) that the pressure is of 
the form p = p + r*-po(¢), so that (0p/dr),-9 = Oand 
po(¢) > 0 along the wall. This means that the inter- 
section of slipstream and wall should occur at a mini- 
mum of the pressure distribution along the wall. This 
feature is confirmed by many of Bleakney’s' interfero- 
grams. 

(C) The vicinity of the singular point S has been in- 
vestigated analytically by resolving the velocities into 
radial and tangential components !/ and V and expand- 
ing the latter in terms of polar coordinates 7, g around 
S. In this way, a recursion formula between the vari- 
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ous coefficients is obtained, which shall be derived in gram suggests. By expansion of U’ and V in powers of A 
the next section. r and substitution, it can be shown that the equations th 
of motion and boundary conditions can be satisfied by le 
AN ANALYTIC INVESTIGATION OF THE FIELD NEAR S_ purely radial isentrops only in the limit of very small fle 
cul , ; . ? ' r—i.e., up to first powers of r—and that for the com (I 
he equations of motion may be written in polar ‘ - 3 i : 
; : : plete solution a more general form of the following W 
coordinates 7, ¢ as follows: ; 
structure has to be assumed :* Ir 
plU,+ (1 ry(V, + U)]+U-p,+(V r)-p, + 2p = 0 ¢1 
(10) e 
a s 7 ra = U = —r*-Cll + > r*-up(¢) | 7 
U-U, + (V/r) (U, — V) + U + (1/p)p, = 0 (11) Pipa" . 
1 
ry. > > / > > > : = 9 = ’ ' 
U-V, + (V/r) (Vy + U) + V + (1/p-r)p, = 0 (12) pe —r*-Cl0 + ¥ r*-x%(¢)] ne 
. : k=l in 
p = (p/o) (13) (14) fo 
ren . ° . . ° a | B u ‘ 
Che behavior of the isentrops near S suggests that, in P = Ptr’: Dr": Prl¢) is 
~ . . > my a y 0 
first approximation, UV = U(r, g) and V = 0 and that, an 
for larger distances " V be different from zero, to ac- _ P | - 3(¢) + r°- 3° 4. p(y) th 
count for the deviation from ray behavior, as the dia- a(r, ¢) n=0 Fe 
ne 
ca 
Substitution of Eq. (14) into Eqs. (10) to (12) yields the following set of equations for the first powers in r: 
Eq. (10); 7°: p(—C — C + 2) = 0, which yields a = 1,C = 1 5, 
r’: p(—3wy — vw’) = 0 or p(—3u — wv’) — w-p'(¢) = O (1) ca 
r?: p(—4u,; — v,’) = 0 or p(—4u, — v,’) — 1 -p’(¢) = O (IT) pe 
. , — ap 
Eq. (11); 7°: no terms in the velocities; hence, 8 + | Re 
> +1—1=6. Since C = 1,86 + 2 ‘i 
r*: 2u + (1/p)-3po(¢) = 0 yields B = 3 (I th 
r3: 3u, + 2m? + r[uo’ — vw] + (1/p)-4p, = O (I] _ 
Eq. (12); 7° andr’: no terms in u%, %, and pp in 
r?: Qu) + (1/p) po’ = O confirms a + 1 = 2 and 8B = 3 (I) 1S 
r*: 37] + SUpVo +- Vo’ +- ( l p) pr’ = Q (II } tre 
alc 
oy . . rin . . . . ‘ . “* . - tic 
Summarizing: The continuity equation yields a = | Solutions: Since on the wall—i.e., for g¢ = ¢, m: 
and C = 1. The first momentum equation yields 8 = —fy'(¢1) = 0 oa 
3; the second momentum equation is also satisfied by F {oe or as 17 qu 
gs > * 7 9 = Aj)SIN og + COtY 5¢g)-COS HE ( 
the indicated values of a, C, and 8. ‘ ” siti =e ' die 
After introducing the dimensionless variables Therefore, inf 
~ . = _ 7 pi\1 l col 
Po = Pole) ai") Pi pi” + 16p; = —18A? ( ) (Is a ¢ 
yo = Up-Qy p/Y sin* 35¢ — 
3 4q 
oe aes with the solution | 
combination of Eqs. (I) yields F 9 -, /pi\l I ee _ 
ply) = —— A? — + C2sin 4g + we 
i i S p/yY sin*® 3¢ | I 
Po (¢) +r 9+ pol) = 0 the 
iy = —(3/2) fo (15) 
i 9 : , coty 4¢)- cos 1p! = 
% = —(1/2)po F ( bot 
: : 7 . ° 4° ~, que 
Combination of Eqs. (II) yields after the boundary condition /,'(¢;) = 0 has been ap kn 
i . rT . . . - . ( 
plied. Thus, the pressure distribution function, near tio, 
}. »F i~ 9 / ad 9 > ~ . O 
Pi” + 16p; + 8[%? + (1/9)%’ 2] = 0 (16) S,is given by 
) 
(r, v) D eS) a eee Pe owas 
t = t + ( ) -A-[sin 38¢ + cotg 3¢,- cos 3¢] + ( ) -[C {sin 4¢ + cotg 4¢,-cos 4¢} — / 
Pi pi ay ay Si 
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Actually, A and C should be determined by purely 
theoretical considerations, from the conditions preva- 
lent on the other end of the flow field—i.e., on the re- 
flected shock. Since the radius of convergence of series 
(19) is CO, the theoretical formula has been compared 
with the experimental pressure distribution on the wall, 


from which A = —2.74, C —2.35; also note that 
¢: = 0.85. 

The isobars = (p — p)/p; = constant of the pres- 
sure field (19) have been plotted in Fig. 2, in which the 
three isobars p = 0 and the alternation of positive and 


negative p-ridges is most conspicuous. This feature is 
intimately connected with the fact that the expansion 
for P starts off with the third power of 7, which in turn 
is a necessary consequence of the equations of motion 
and the “‘ray behavior’ of the isentrops. Except for 
the three isobars p = 0, no other isobar can originate in S. 
For finite values of /, ry can never become zero because 
neither of the two trigonometric expressions in Eqs. (19) 
can grow beyond bound. 

The behavior of the isopyenics p/p) = constant, near 
S, is in striking contrast to that of the isobars. Be- 
cause of Eq. (14), p assumes a different value in S, de- 
pending on the direction ¢ from which this point is 
approached; in other words, the zsopycnics demonstrate 
a similar “ray behavior’ as the isentrops. Taking ac- 
count of the fact that the range of o-variation must be 
the same near S as on the reflected shock, the range of 
isopyenics emanating from this point can be computed 
in accordance with Eq. (14). Experimental evidence 
is scarce, since the p- and p-field in this region is ex- 
tremely flat. However, since the main change of o 
along curve NT is concentrated at the ends,' in par- 
ticular near 7’, it must be expected that at S relatively 
many isopyenics arise near wall and slipstream and 
only a few in the center sector, as indicated in the 
qualitative isopyenics diagram Fig. 3. At the occa- 
sion of a personal, visual evaluation of Bleakney’s 
interferograms, the authors believe to have found a 
confirmation of Fig. 3, which could be improved into 
a quantitative diagram, similarly as Fig. 2, by use of 
Eqs. (10) to (13). 

It may have been noted that in the preceding treat- 
ment of the field near S, a boundary condition was 
prescribed only on the wall, but not on the slipstream; 
the latter may be identified with any one of the isen- 
trops, so that its shape is determined by the other 
boundary condition and the additional requirement of 
quasi-ray behavior. Since the slipstream form is not 
known a priort but should be determined by the solu- 
tion itself, this result is quite satisfactory. 


COMPARISON OF THE FIELDS NEAR S AND O 


A further elucidation may be obtained, if the field near 
Sis compared with the field in the vicinity of the wedge 
corner O, where another type of singularity arises, 
namely a true stagnation point. Also here Eqs. (10) 
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to (13) apply, but now the requirement of tangential 
flow direction has to be enforced on both walls; on the 
other hand, no ray behavior of the isentrops may be ex- 
pected near O, as clearly indicated by the diagram. 
Therefore, it is to be required that T(x) = T(38°) = 0, 
and that, in contrast to Eqs. (14), ( and VP start off 
with the same power of r. Furthermore, near the stag- 
nation point, the field must go over into the incompres- 
sible solution which satisfies the equations of motion, 
without the terms U’-p, + (VT r)-p, in Eq. (10) and 
with p = const. instead of Eq. (13). This yields 
around 0 the solution : 


U = A-r*-. cos [k(r — 6)] — 1 
v= —A-r*—. sin [k(9r — 0)] > (20) 
p=p- r?®*—2.(p4/2) + 
Both boundary conditions (no circumferential velocity 
on either wall) are fulfilled, if k = 1.28. Eqs. (20) 


gives a unique starting point for the series expansion 
represeiiting the compressible solution, which can be 
derived by use of the complete equations of motion. 
It is seen that, in the limiting incompressible case, p = 
const. for 7 = const., 1.e., the zsobars are circles; and that 
in the next approximation, p = const. for r = const., 
i.e., the tsopycnics are circles around O, as confirmed 
by the experimental diagram, and in contrast to the 
isopycnics near 5S. 

The preceding analytical discussion is in agreement 
with the results obtained by numerical integration, 

Summarizing the comparison of the fields near S and 
O, one may state: In either case, Eqs. (10) to (13) 
have to be fulfilled, and either solution has to go over 
into the respective incompressible solution near the 
origins. Both boundary conditions can be fulfilled 
only for fractional, leading power of r (i.e., fractional 
k). On the other hand, if & is fractional, the equations 
of motion cannot be fulfilled, unless both velocity 
components start with the same power of r; in other 
words: ‘“‘Ray behavior’ of the isentrops is incom- 
patible with the fulfillment of both boundary condi- 


tions. 


THE ALTERNATIVE, HYPERBOLIC METHOD 


The theoretical difficulties which we had to overcome 
in the preceding solution procedure, stem, to a great 
extent, from two sources: 

(a) From the nonlinearity of the four differential 
equations of motion (see Eqs. (1) to (4)), which in a 
difference scheme lead to a large system of equations 
of the second and third degree, and 

(b) From the “floating character’ of the boundary 
conditions on the shocks and slipstream which, in our 
numerical scheme, lead to a trial and error procedure, 
until the correct shock location is found. 

It is to be noted that NV equations of the third degree 
have 3” sets of solutions, and care must be taken that 
any iteration procedure converges toward the physically 
correct set. This as well as the labor involved in de- 








termining the shock location makes some experimental 
advance knowledge of the results highly desirable in 
the preceding method. 

On the other hand, since it is often required that 
the theory predict new phenomena, it is a worth-while 
objective to develop, if possible, a self-consistent theory 
which operates without the aid of experimental hints. 
An idea promulgated by J. von Neumann? and carried 
further recently by P. Lax® suggests a scheme in which 
both difficulties mentioned under (a) and (b) can be 
avoided. 

As a starting point, the equations of motion are not 
expressed in conical coordinates, as is done in Eqs. (1) 
to (4), but are left in their original x, y, ¢-coordinates, 
in which their character is always hyperbolic, no matter 
what kind of flow field is considered; see Eqs. (21) to 
(24). This is essential, because it is this feature which 
allows the floating boundary conditions to be removed. 
Applied to the problem of Fig. 1, it is completely suffi- 
cient to impose one boundary condition only on the 
wall. If the initial field at ¢ = 0 is prescribed, the solu- 
tion everywhere in the whole x-y-plane is determined 
for all later times, and the shocks and the slipstream 
will come out automatically, in the form of relatively 
sudden steep drops of density, pressure, etc., at the 
theoretically correct location. In other words: The 
numerical step-by-step procedure, which is adapted 
to the hyperbolic character of the equations, yields the 
shocks as part of the general solution defined in the 
whole infinite plane. 

While in this manner the difficulty mentioned under 
(b) is removed, the nonlinearity indicated under (a) is 
rendered harmless at the same time: One can see that 
in the equations of motion, if written in hyperbolic 
form (see Eqs. (21) to (24)) all time derivatives appear 
linearly, while the space derivatives occur in nonlinear 
combinations. Because of the preferred status of the 
time coordinate, ¢-derivatives shall be represented as 


forward differences :° 


: I w 
fy = J pnt sd 


At 
Nn ‘au rit pt 5 
fi + 1, m 7 Fi 1, m nei? m+ 1 aE m if (A) 
4 f 
space derivatives, however, by central differences: 
n cn 
- Ii + l,m ee 1, m (B) 


2-Ax 


here superscripts refer to the time cycle, subscripts to 
the position in space. After substitution of these 
expressions, it can be seen that at the instant » + 1, 
the only unknown f;' ** appears linearly, while any 
nonlinear combination of the other quantities f; are 
taken at the previous instant m, at which conditions are 
supposed to be known. 

Taking into account the preceding remarks, we have 
formulated below the analytical statement of the shock 
diffraction problem represented in Fig. 1: 


Equations of Motion’ 


pic t+ div M=0 (21 
m,+ div R; = 0 (22 
r, + div Rk, = 0 (23 
E, + div T = 0 (24 


here p is the density, MW the mass flow vector with the 


components m and r, and £ the total energy (sum of 
kinetic and internal); furthermore R,, 7» 
vectors, which are certain 


E 


R, = i Gj- E+ "— (") - 
2 \p 


M7 and p, namely: 


R, = [| + lo = je 4 (25 
p 


T = (= ree Bees yim a. 
p 9 p” 


Now, write, e.g., Eq. (21) in difference form: 


1 At | F 
— mM) 
2 Ax 


1 At | ? P 
r; — Fy 3 1 (26) 
2 Ay 


Since all, right-hand terms are taken at cycle ¢, they may 
be considered known, and at time (¢ + 1), pid ' can be 
computed. Of course, At/Ax and At/Ay are related 
to the characteristics of the hyperbolic system Eqs. 
(21) to (24), and have to be chosen in accordance with 
the stability criterion of Courant, Friedrichs, and Lewy.’ 


Boundary Condition 


ris, ¥,. $) Vv . = 
: = - = given tang 86, on wall (27) 
m(x, y, t) E 
Initial Condition 
For z <0, 7 > 0: ) 
m(x, y, 0) = given const.; r = 0 | 
p(x, vy, 0) = pi; E(x, vy, 0) = E, 
For x > 0, y > 0: (28) 
m(x, ¥, 0) = 0; r=0 
. 0 
p(x, vy, 0) po; E(x, y,0) = - } 
y¥-—! 


and F are 
nonlinear combinations of 
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As seen from the initial condition Eqs. (28), the flow 
ficld at ¢ = Ois that of a plane, undisturbed shock front 
approaching a wedge. 

After writing all four Eqs. (21) to (24) in finite dif- 
ference form, it must be expected that by evaluation of 
these equations at successive time instants, the whole 
diffracted shock field, including slipstream and shock 
location, will gradually emerge through successive iter- 
ations, in a movie-like fashion. 


JUSTIFICATION OF THE HYPERBOLIC METHOD AND FIRST 


PRELIMINARY RESULTS 

On first sight, the above expressed expectation may 
seem too optimistic to the reader who is not familiar 
with the character of hyperbolic equations and differ- 
ence schemes. And, indeed, it is most surprising that 
a so simple configuration, as a plane incident shock 
field, can without extraneous guidance, be transformed, 
by step-by-step computations, into a so complicated 
shock field as that of a Mach configuration. How- 
ever, it must be kept in mind that, if some equivalent 
of an existence and uniqueness theorem exists,® it must 
be possible to construct the solution of the initial- 
boundary-value problem by systematic application of 
the correct equations of motion, boundary condition, 
and initial field. 

Now, conditions in this problem are somewhat com- 
plicated by the fact that the initial field contains a dis- 
continuity surface (shock), across which the equations 
of motion (21) to (24) do not hold, but have to be re- 
placed by the conservation laws, in which an entropy 
difference appears. On the other hand, it has to be 
realized that in a difference scheme, the two points on 
both sides of the same shock surface are apart at least 
by one (or several) mesh distances, and that therefore 
the distinction between equations of motion and con- 
servation laws becomes meaningless; in other words: 
In difference form, the equations of motion are equivalent 
to both relations, and hold, without exception, throughout 
the whole field of flow. 

Thus in the gradual build-up of the Mach field, two 
phases may be distinguished: first, the transforma- 
tion of the discontinuous initial field at ¢ = 0 to an en- 
tirely continuous field (in difference form) after a very 
few cycles. During this first phase, a geometrically 
dissimilar development takes place, namely the start 
of the formation of a compression wave at the wedge 
corner O. Second, there follows the transformation 
of this latter continuous field into the final Mach con- 
figuration, after very many cycles. This involves a 
development which might be called asymptotically 
similar. 

The application of the difference formalism erases 
automatically the initial shock discontinuity and pro- 
duces the desired continuous field, after the first few 
From there on, a continuous initial value prob- 
The fulfillment of the one 


cycles. 


lem has to be solved. 
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boundary condition on the inclined wall produces auto- 
matically a (nearly circular) reflected compression 
wave, similarly as the boundary condition on a plane 
wall produces a plane reflected wave. In subsequent 
cycles, the slope of the reflected wave will gradually 
steepen, similarly as a plane pulse is known to steepen 
and finally form a shock. In general, this phenomenon 
is treated as an isentropic process, while in our case a 
considerable entropy difference must develop across 
the shock. This is taken care of by the application of 
the energy equation. Again, in analytic form, the en- 
ergy equation is equivalent to DS Dt = 0; and any 
particle keeps its entropy along its streamline. But, 
in difference form, the energy equation is equivalent 
to both, the equation of motion as well as the corre- 
sponding conservation law, and therefore it is per- 
missible to impose, initially, a given steep entropy 
and let it 


spread and then steepen up, by application of the step 


gradient (i.e., across the incident shock) 


by-step procedure. 

Moreover, the following interesting feature is to be 
noted. Although the motion of an ideal fluid was as- 
sumed, gradually S-values outside of the initial range of 
S-values appear, although no account is taken of any 
viscosity term which could be the source of such new 
entropies. According to a suggestion by K. O. Fried- 
richs, this may be ascribed to the difference in form 
between the differential and difference equations of 
motion ; the difference between corresponding terms may 
be interpreted as artificial viscosity terms. 

More explicitly: If each difference quotient in the 
equations (which is of first order) is expanded in powers 
of a mesh width parameter, the set of first terms in the 
expansion represents the differential equations of an 
ideal fluid, and the set of second expansion terms (which 
are second order derivatives) represents the artificial 
viscosity terms (similarly as in the Navier Stokes equa 
tions). Although the amount of viscosity is not quite 
the true one, the final result will be correct, since the 
final situation may be described in terms of ideal fluid 
concepts only, as the Rankine Hugoniot equations 
show. And any viscosity term will merely act as “‘cata- 
lyzer’’ toward that end. 

From the preceding it can be seen that the final cor- 
rect field can be obtained only after many cycles, 1.e., by 
the use of modern electronic computers. We have 
coded the described problem for the I.B.M. ‘Defense 
Calculator,” for a shock strength 10:1 and a wedge 
angle of 60°. However, in an attempt, by disregarding 
any details, to merely verify the essential features of 
the build-up of the Mach configuration, we have carried 
out the first six cycles by hand computation. The re- 
sults are shown in Figs. 4 and 5. 

Fig. 4 represents the density distribution along the 
wall and symmetry line, during the first six cycles. 
One will note that during that time interval a broaden- 
ing compression wave forms at the wedge corner, and 


that the steep density gradient, representing the 
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DENSITY DISTRIBUTION ALONG WALL IN SUCCESSIVE TIME CYCLES 
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Mach stem, moves up the inclined wall, while the left 
slope of the compression wave moves very little. This 
is in accordance with the fact that the compression is 
continuously swept to the right, in the wake of the 
incident shock. Furthermore, the stagnation density, 
which alternates about an equilibrium value, is of the 
correct order. 

In Fig. 5, the density distribution in polar view is 
shown. The shaded areas are the regions of maximum 
density gradients. One may see that already in the 
4th and 5th cycle the correct split-up of the Incident 
into Reflected and Mach shocks has occurred. 

Naturally, details are still missing here, and may be 
expected only after an order of about 50 or 100 cycles. 

Note: It has been pointed out that the develop- 
ment of the final flow field consists of two phases: a 
first phase during which geometrically dissimilar 
changes take place, and a second phase, during which 
the changes are asymptotically similar. The preceding 
hyperbolic method seems to be adequate for treating 
the first phase (during which a new configuration 
forms). But in the second phase a procedure em- 
phasizing the conical character of the emerging field 
would be preferable. Now, both requirements can be 
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met, by transforming the equations of motion from an 
x-, y-, t-system into a é-, n-, t-system (here & = x/t; 
n = y/t), and separating the ¢-dependent terms from the 


remainder: 
div M+ E-p: + np, = t-(pide y 
div R; + E-m; + n-m, = t-(mi)z , 


div R, + &-r; = n°T, = t- (ride n 


div 7+ &-E,+7-E, = t-(Ed;, 


In this new form, nothing is lost of the generality of 
the equations, but an opportunity is afforded for em 
phasizing the conical character, if the problem is of this 
nature. 

If the time derivatives on the right side are disre- 
garded, the equations degenerate into those of a conical 
field, and the representation of the flow field becomes 
static. But so long as the time derivatives are differ- 
ent from zero, the problem is hyperbolic in time, and 
the step-by-step procedure can be applied as before. 
Hence, the values of all time derivatives must become 
automatically smaller and smaller during the emergence 
of our final field. Their magnitude is a criterion of 
how far we are away from the final goal. Now, it ap- 
pears possible to accelerate the convergence of this 
process by giving the conical terms in the differential 
equation a heavy weight, after the first few cycles. 
The introduction of such a weight factor is tantamount 
to adjusting the expressions Af/Aé and At/An in such a 
fashion that the time derivatives become small after as 
few steps as possible. How this can be arranged with- 
out violating the stability criterion is at present under 


study. 
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A Simple Laminar Boundary Layer 


with Secondary Flow 


* 


HENK G. LOOS+ 
California Institute of Technology 


SUMMARY 


The incompressible laminar boundary layer over a flat plate 


is studied for the simple case where the stream lines in the free 


flow have a parabolic shape. 


layer equations is derived. No separation occurs, even when 


An exact solution of the boundary- 


there is a strong adverse pressure gradient along the streain lines, 


so that in this instance the secondary flow has a favorable in 
fluence. 
stream line to another in the free stream, the total pressure within 
the boundary layer at a given point can exceed that of the cor- 


responding free stream 


SYMBOLS 


hs = Cartesian coordinates 

d u = velocity components in x, y, 2 directions, respec 
tively 

U,V, W = velocity components u, v, w in the free stream 

1 = magnitude of the total velocity in the free stream 

tn = velocity components respectively tangential and 
normal to the free flow direction 

6 = angle between free flow direction and the normal 
to the leading edge in an arbitrary point 

p = static pressure 

v = kinematic viscosity 

a,b = constants 


INTRODUCTION 


Ws AN AIR STREAM FLOWS over the surface of a 
flat plate, the curvature of the stream lines 
within the boundary layer may differ from that of the 
free stream because of a gradient of pressure normal to 
the direction of free-stream flow. Then, since accord- 
ing to boundary-layer theory the pressure in the boun- 
ary layer does not vary with distance from the surface, 
the boundary-layer flow must, because of its lower ve- 
locity, curve more sharply than the free-stream flow in 
order to balance this pressure gradient. The boundary- 
layer velocity component normal to the direction of the 
free stream has become known as the cross flow or the 
secondary flow. This secondary flow arising from a 
boundary layer operating in an externally produced 
pressure gradient was initially discussed by Prandtl! 
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Because of the variation of total pressure from one 
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and is of considerable importance to the understanding 
of complex flow phenomena that occur in turboma- 
chinery. 

A detailed investigation of the cross flow in the 
laminar boundary layer on an infinite swept wing was 
carried out by Sears,” where the pressure gradient in 
the free stream was normal to the leading edge and, 
hence, had a component normal to the free-stream 
flow. Sears calculated the resulting cross flow and 
showed in particular that the cross flow vanishes for the 
swept flat plate at zero angle of attack that is, when 
the imposed pressure gradient vanished. 
flow generated by the centrifugal pressure gradient on 


The cross 


a propeller or turbomachine blade was investigated in a 
similar manner by Fogarty,’ where the blade was 
treated as an airfoil of infinite span rotating about a 
point of its span. 

The study of three-dimensional boundary layers has 
been pursued further by Howarth‘ and other workers. 
A problem of particular interest in the study of turbo- 
machines was discussed by Mager and Hansen,’ who 
considered the cross flow generated by a curved free 
stream, carrying its own pressure gradient, which flows 
over a semi-infinite flat plate. When the free stream 
is a potential vortex, a solution was obtained which is 
valid over distances from the leading edge for which the 
turning angle of the free stream is small. 

Now it seems clear that a considerable portion of the 
cross flow in boundary layers is directly related to the 
fact that the velocity in the boundary layer is less than 
in the free stream. For example, if a thick boundary 
layer is built up by a straight uniform flow over a plate 
and then at some point the free stream is turned, the 
cross flow that takes place has little to do with local 
viscous stresses but is simply a matter of transporting 
the vorticity that was generated far upstream. This 
idea was used by Squire and Winter® to discuss the 
secondary flow that occurs in airfoil cascades having a 
thick wall boundary layer and by Hawthorne’ in his 
treatment of the flow in pipe bends. The process 
becomes simpler when viscous stresses may be neg- 
lected, and it has often been suggested (e.g., Hayes*) 
that the general problem of three-dimensional boun- 
dary-layer flow may be considered in two parts: an in- 
viscid outer part matched with a viscous sublayer. 

The present work, however, is an exact calculation 
of the laminar boundary layer developed on a flat plate 
by a flow having stream lines parallel to the plane of the 
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plate and of parabolic shape in this plane. The free 
stream is rotational, having a constant vorticity that is 
directed normal to the plate. This situation is similar 
to that which occurs in some problems of fluid motion 
past turbomachine blades. 

MAIN STREAM 


BOUNDARY LAYER WITH PARABOLIC 


Let U, V, and IV denote the components of free- 
stream velocity in the x, y, and g directions, respec- 
tively, where x is the distance measured along a semi- 
infinite flat plate in a direction normal to its leading 
edge, y normal to the plane of the plate, and z normal 


to these two coordinates as shown in Fig. |. If the 
velocity components are chosen to be 
l” = constant 
Yo =f (1) 
We = il + by 


with a and 6 constants, then the stream lines of the 
main flow are of parabolic shape and the flow is of uni- 
form vorticity, 2 = 6, with the vorticity vector in di- 
rection of the y axis. The Helmholtz relation for for- 
ticitv is obviously satisfied. 

The boundary-layer equations for an incompressible 
fluid are, in this case, where the velocity is independent 


of g, 
uO 4 Ou Oru (9 
u v = yp 2) 
Ov ov Ov" 

Op oy = 0 (73) 

Ow Ow : -w 
u - v —bU=py = (4) 

On’ Ov Ov" 
(Ou Ox) + (Ov/O¥) = 0 (5) 


with the boundary conditions 


u=0 
y = 0, v =0 (6) 
w= 0 
v= 6, = ay = 
i \ (4) 


lw=at bx 


Eys. (2) and (5), together with the boundary conditions 
for the velocity u, are identical with those for the two- 
dimensional laminar boundary layer on a flat plate with 
Therefore the solution for 1 


zero pressure gradient. 
Eq. (4) and 


and v is the well-known Blasius solution.’ 
the boundary conditions for w remain for determination 
of the velocity component w. 

Eq. (4) can be reduced to an ordinary differential 
equation by using » = y+/U/vx as new independent 
variable. The Flasius solution may be expressed in 
the form 


“w= TF’) (S) 


Vx 
v= - - g(n) (9) 
x 
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my 
g(n) = | mi" (m) do (10 
Therefore if the velocity component parallel to the 
leading edge is assumed to be of the form 
w= Wy (7) + bxw,(n) (11) 
Eq. (4) may be transformed to 


—(1/2) (F'n — g)iey’ — wo"] + 
b — | + F’w, — (1 2) x 


( fr’ so 


/ ” . 
2 g—~ 27a, ~~ BW = 0 (12 


(a/x) 


Furthermore, by partial integration of Eq. (10), the 
quantity F’y — g appearing in Eq. (12) may be written 


F'n ss F'(m) dn 


so that Eq. (12) becomes 
(a/x) [—(1/2) Fao’ — wo"”] + 
b{ —1 + F’w, — (1 ?) Fwy,’ one Ww,” = 0 (13) 
Now if the form assumed in Eq. (11) is valid, Eq. 
(143) must be an identity in v and the quantities wy and 
w fulfill the equations 
wy” + (1/2)Fu’ = 0 (14) 
Ww,” + (1 2) Fwy,’ = F'w, = —] (1.5) 


with boundary conditions 


= U0, wy = O 
0 ‘ (16) 
ww) = 0 
= adi wy = l al 
” (17) 
Ww) = l 


The Blasius solution for the velocity u = UF’(n) in the 
direction of x satisfies the differential equation F”F + 
2F’”’ = 0, with the boundary conditions 7 = 0, F’(0 
n = ~, F’ = 1. It is obvious that the function wo = 
F’(n) satisfies Eq. (14) and the boundary conditions 
[Eqs. (16) and (17)]. The cross-flow component wp, 
which is related with the constant free-flow component 
IV = a, has the same profile as the velocity component 
u. This result was first obtained by Sears in reference 
2. It means that in the case of a parallel stream ap- 
proaching a flat plate at an arbitrary angle, the flow 
direction in the boundary layer is the same as in the 
free stream. In other words, no cross flow occurs, 
and, consequently, a cross flow in a boundary laver can 
only exist if there is a pressure gradient. 


CALCULATION OF THE CROSS FLOW 


The Eq. (15) is linear and nonhomogenous; the 
coefficients F(y) and F’(n) are tabulated functions.° 
Instead of calculating the solution of Eq. (15) with 
boundary conditions given by Eqs. (16) and (17), « and 
w will be related with velocity components ¢ and n, 
which are, respectively, tangential and normal to the 
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local free-flow direction. In Fig. 2 a free-flow stream 
line is shown for the case that the coefficient a is zero 
that is, so that the stream lines are perpendicular to 
the leading edge of the plate. If # is the angle between 
the free-flow direction and the positive x axis, 


t= ucosd8+wsind (18) 

n = weosd — usind (19) 

Now, for a = 0, u = UF'(n) and w = bxw,(n), and 
/ 

T = VU? + 6*x? is the local free-stream velocity. 


Then Eqs. (18) and (19) may be written as 


t uF’ + b?x?w, F’ + (bx/u)?w, 


= , = oe (20) 
7 u? + b*x? 1+ (bx/U)? 
n bx Uw, — bx UF’ bx w, — F’ (21 
= = 21) 
i U2 + bx? U 1 + (bx/U)? 
Now it is convenient to define 
h(n) = wi(n) — F’(n) (22) 


Then a differential equation for h(n) can be derived in 


the following way: The function w;(7) is a solution of 
w,” + (1/2)Fw,’ — F’w, = —1, and F(n) satisfies the 


equation F’’’ + (1/2)FF" = 0. Subtraction of these 


differential equations gives 
(w, — F’)” + F(w, — F’)’ — F’(m, — F’) = 
—1 + (F’)? 


and taking account of the definition of h(n), 
h"(n) + (1/2) Fh’ — Fh = —-1+ (F’)? (23) 

The function h(n) is equal to —G’(n) which has been 
tabulated by Mager and Hansen.’ The normal or 
cross-flow component |Eq. (21)], may now be written 
as 

n bx h(n) 

= = (24) 

T U 1+ (6x/U)? 
For small values of xv, (6x/U)? <1, and the normal 
flow component is the same as that found by Mager 
and Hansen’ for small turning angles. This result is 
obvious, because, for a flow with a = 0, the circular 
stream line coincides in first order with a parabolic 
stream line for points close to the leading edge. 

From Eq. (24) it follows that for large turning 
angles in “‘parabolic flow” the distribution of the rela- 
tive normal flow component 2/7 remains the same and 
the magnitude varies as (bv/U)/[1 + (bx/U)?]. It is 
possible to express n/7’ in terms of the turning angle 
v? of the free flow. For a = 0, bx/U = tan @ and 


(bx/U)/[1 + (bx/U)?] = (1/2) sin 28 
and therefore 
n/T = (1/2) h(n) sin 28 (25) 
Further, Eq. (20) can be written as 


t/T = F'(n) + h(n) sin? 3 (26) 


In the general case where a ~ 0, w = awn) 4 
bw;(n), the normal velocity component can be expressed 


as 
n = [aF’ + bx(F’ + h)| cos 3d — UF' sin d 
Now, tan 3 = (a + bx)/U, so that 


n = bxh(n) cos 3 (27 


Further, the quantity bx may be written bx = U(tan 
dv — tan J), where & is the angle between the free-flow 
direction and the x axis at the leading edge (see 
Fig. 3). Eq. (27) becomes 

n/T = A(d, do)h(n) (28 
A(d, Jo) = cos? d(tan F — tan Vo) (29 


In the same way an expression for t/7’ may be derived 
Ald 


in terms of F’(n), h(n), & and J, 
t/T = F'(n) + BV, do)h(n) (30 


B(d, Jo) = (1/2) sin 20(tan 3 — tan Vo) (31 
RESULTS 


The functions A(#, Jo) and Bd, Jo) are plotted in 
Figs. 6 and 7 for d) = —60°, 0°, and 60°. For the same 
values of J, the tangential flow ¢/T is plotted versus 
n for some values of Jo. 

In Fig. 8, % = —60° and 3d = —60°, —40°... 60°. 
For 0 = —40° and 8 = —20°, 
profile is found, caused by the adverse pressure gradient 
along the direction of the free-stream flow. As will 
be shown later, the cross flow prevents any actual sepa 


a ‘‘separation”’ type 


ration. 
After passing the point 3 = 0, the velocity profile 
becomes fuller and at the point 3} = 20° the velocity 


in the upper part of the boundary layer exceeds the 
free-stream velocity by a small amount. For greater 
values of J, the maximum velocity in the boundary 
layer increases more and more, until at 3} = 60° the ve 
locity distribution with the highest velocity peak is 
reached. For 3 > 60°, the maximum velocity in the 
boundary layer drops again, and in the limit 3) = 90 
the velocity distribution is just F’(n) + h(n), still hav- 
ing an increased velocity in the upper part of the boun- 
dary layer. For d& = 0° and & = 60° (Figs. 9 and 10) 
the velocity profile becomes fuller with increasing 
dv, and the highest velocity peak in the boundary layer 
is reached at 3} = 90°. Of course, for this value of , 
no “‘separation”’ type profile occurs because there is no 
adverse pressure gradient for the tangential flow. The 
normal or cross flow n/T is distributed over n as A(d, 
Jo) h(n), and so the velocity profile will be always of the 
shape /i(n) (see Fig. 5) while the amplitude changes as 
A(#, Jo) (see Fig. 6). 

It is of interest to investigate the total head distribu- 
tion in the boundary layer. Because, in the boundary- 
layer approximation, the pressure is constant through 
the boundary layer, the total head loss in the boundary 
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layer can be determined by calculating the kinetic 


energy. Neglecting the velocity component v, the 


ratio of the kinetic energy in the boundary layer and in 


the free flow is, for 3) = 0, 

(n) (UF’)? + (F’ + h)?b2x? 

E( ) U? + bx? 
and with 
U/V U2 + bx? = cos J, bx/V U2 + bx? = sin J 

E(n)/E(~) = (F’)? cos? 38 + (F’ + h)? sin? 3 

If the ratio /{(n) /:(@) exceeds unity, the total head 
at the point of the boundary layer exceeds the free- 
stream value. If 7 = 2.0 is chosen, according to Fig. 


5, F’(2.0) = 0.63 and F’(2.0) + h (2.0) = 1.17. Thusif 


E(2.0)/E( 2) = 0.63? cos? 8 + 1.177 sin? 8 > | 


that is, if } > 52.1°, the energy in the boundary layer 
at » = 2.0 exceeds the free-stream energy. 

This result can be understood by comparing the 
stream lines of the free stream and of the boundary 
layer (Fig. 4). Since the stream lines in the boundary 
layer curve more sharply than those in the free stream, 
the stream line passing the observation point P at 7 = 
2.0 (dashed line in Fig. 4) comes from another upstream 
location than the free-flow stream line passing through 
the same point P. There is a total head gradient in 
the z direction upstream of the leading edge because of 
the vorticity in the parabolic flow, so that upstream 
of the leading edge the dashed stream line through P 
will have a larger total head. Along the dashed stream 
line, energy will be dissipated, but it is possible that the 
energy loss is less than the difference in total head up- 
stream of the leading edge. In this case, the total head 
in some point of P will exceed the free-stream value. 

Another interesting question is whether there will be 
a separation of the boundary layer. Of course, sepa- 
ration cannot be expected for J) > 0 because the pressure 
gradient in this case is favorable. But if —vp is sufli- 
ciently large, an adverse pressure gradient occurs in the 
region 3) < 3 < 0, and if there were no cross flow, the 
boundary layer would separate. Actually, a cross 
flow occurs, and the separation may be studied by 
employing the solution for the boundary-layer flow 
which has been obtained. 

If three-dimensional separation occurs, then, ac- 
cording to Hayes,’ there is a separation line on the sur- 
face which is tangent to all ‘‘stream lines’’ at the sur- 
face. In the present example, the velocity field is con- 
stant in the z direction, so that any separation line must 
be parallel with the z axis. The flow direction 3’ on 


the surface is given by 
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al’ + bx(F’ +h 


tan 3’ = lim = lim a 
yoo U y—>0 UF 
a+ bx [1 + lim h(n) /F'(y 
10 
7 


Now 


lim h(n) /F’(n 

70 
is finite so that, for finite values of x, tan &’ will be finite. 
Therefore no separation occurs for finite values of 
From this result it follows that the cross flow avoids 
separation in the present example by transporting fluid 


with high energy to the ‘‘critical’’ regions. 


CONCLUSIONS 


The boundary-layer equations for the “‘parabolic’’ 
flow along a flat plate can be solved exactly, giving in- 
formation about the secondary flow (cross flow). The 
solution for dj) = O coincides with that of Mager and 
Hansen’ for small turning angles but also remains valid 
for large angles. There, the formulas of Mager and 
Hansen can be used, replacing the turning angle 3 by 
(1/2) sin 20. 

If —d is large enough, a “‘separation”’ type profile for 
the tangential flow will be found, but there will be no 
separation because the secondary flow transports fluid 
with high energy to the “‘critical’’ regions. According 
to the total head gradient upstream of the leading edge 
(rotation), the total head in the boundary layer can 


exceed the free-stream value. 
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Some Considerations on the Air Forces on a 
Wing Oscillating Between Two Walls for 


Subsonic Compressible Flow 


DONALD S. WOOLSTON* ann HARRY L. RUNYAN* 
Langley Aeronautical Laboratory, NACA 


SUMMARY 


rhe problem of the determination of the air forces on an oscil- 
lating airfoil between plane walls has, until recently, been treated 


only for incompressible flow. The present paper is concerned 


with the important effects of compressibility, which may be of 


significance in such problems as the measurement of oscillating 
air forces or of wing flutter characteristics in wind tunnels, and 
in the flutter of airfoils in cascade. The possibility of the exist 
ence of an acoustic resonance phenomenon under certain critical 
conditions is discussed. The integral equation for the com- 
pressible case, as obtained by Runyan and Watkins in NACA 
TN 25852, is reviewed briefly and a method of solving the equa- 
tion is given. The procedure is applied to a number of selected 
cases at various Mach numbers and tunnel heights. The effect 
of the presence of the walls is shown to be very significant near 
the resonant frequency and, for certain conditions, to be large 


even at frequencies well removed from resonance. 


SYMBOLS 


= velocity of sound, ft. per sec. 
1, = coefficients in series expression for 
lift distribution (Eq. IT) 
= wing half-chord, ft. 


h = displacment of wing in vertical 
translation 
I = nondimensional tunnel height 


based on half-chord, } 
H,,? = Hankel function of the second kind 
k = reduced frequency, bw/U 


A(.M, sz) + ACM, 2, H kernel of integral equation 

L(x), Lo) = lift distribution, lb. per ft. per unit 
span 

E. = magnitude of the lift for wing in a 
tunnel 

B = magnitude of the lift for wing in 
free air 

V = Mach Number, l’/a 

p = MWkiI/2rB 

? = ratio of tunnel width to tunnel 
height 

R, = V(x — Xo)? + BA(ni7)? 

l = stream velocity in chordwise di 
rection, ft. per sec 

\ = vertical induced velocity (per- 

pendicular to chord), ft. per sec. 

\ = axis of rotation measured from 
midchord, positive rearward, 
based on half-chord 

C; 2 = nondimensional chordwise  coor- 


dinates, based on half-chord 


Presented at the Aerodynamics Session, Twenty-Second 


Annual Meeting, IAS, New York, January 25-29, 1954. 
* Aeronautical Research Scientist, Dynamic Loads Division. 


y = nondimensional vertical coordi 
nate, based on half-chord 

7 = k(x — x 

a angular displacement of wing in 
pitch 

8 = VY1—- Vf 

€ = ix — xo'/BH 

k = mass ratio (rpb?/m, where wm is 


mass of wing per unit length 
u = \Wk/p? 
0 = —cos vr, 6 O08 a 
= circular frequency of oscillation, 


rad. per sec 


| 
| 


w,, = circular frequency at resonance, 
rad. per sec 

= natural circular frequency of wing 
in vertical translation 


= natural circular frequency of wing 


Wa 
in pitch 

p fluid density, slug-ft.* 

o = phase angle between lift force and 
position for wing in a tunnel 

rr = phase angle between lift force and 
position for wing in free air 

Subscripts: 
n = (Q, 1, 2, 


INTRODUCTION 


Shes INFLUENCE OF constrained flow, existing, for 
example, in a wind tunnel, on the forces acting on 
an airfoil has been a continuing problem for many years. 
In the case of steady flow the subject has been exten 
sively studied and reported. In general, the investiga 
tions of the steady case have led to the determination 
of relatively simple factors which can be used to modify 
measurements of the air forces on a wing in a tunnel to 
correspond to free-air conditions. The extension of the 
steady flow results to compressible flow conditions 
presents no serious difficulties (except at transonic 
speeds) since the results for incompressible flow can 
be modified according to Prandtl-Glauert type correc- 
tion factors. 

In contrast, however, the corresponding problem of 
the effect of walls on an oscillating airfoil has received 
relatively little attention, particularly in the case of 
compressible flow. Iniicompressible flow, the problem 
has been studied by W. P. Jones,' E. Reissner,’ and R. 
Timman.* In these papers, the influence of the walls 
is found to be relatively small, although indications 
are given that for some combinations of values of re- 
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duced frequency parameter and a certain tunnel geom- 
etry the effect may be significant. The effect of walls 
in compressible flow has been treated, for the two- 
dimensional case, by Runyan and Watkins,‘ who dem- 
onstrated that compressibility, in addition to its role 
of modifying the forces in the usual sense can also, 
under certain conditions, result in a transverse reso- 
nance phenomenon which can greatly change the nature 
of the flow. This arises from the fact that, in an in- 
compressible fluid, the velocity of propagation of a dis- 
turbance is infinite and no time lag occurs between the 
initiation of a disturbance and its effect at any other 
position in the field, but, in a compressible fluid, a 
definite time is required for a signal to reach a distant 
field point so that both a phase lag and a change in 
magnitude result. Under certain conditions this phase 
lag can give rise to a resonant condition which would 
involve large corrections. 

In the paper of Runyan and Watkins‘ an integral 
equation which relates the motion of a wing to the load- 
ing on the wing was formulated for the case of an os- 
cillating airfoil between walls. In this earlier work, the 
condition of resonance was shown, but other significant 
details of the problem were not delineated. The pur- 
pose of the present paper is to continue the investiga- 
tion of the effects of compressibility by indicating a 
procedure for solving the integral equation and by 
applying this procedure to a number of selected ex- 
amples. The examples treated show the effects of 
tunnel walls associated with variations in tunnel speed 
or Mach Number, frequency of oscillation of the wing, 
and tunnel height. In addition, a brief investigation 
is made of the problem of flutter of airfoils in cascade 
for the case of potential flow. 

The calculation procedure and the results contained 
in this paper are of significance for such problems as 
the experimental measurement of the forces on an air- 
foil and the determination of wing flutter character- 
istics in wind tunnels, and also in certain possible types 


of flutter of airfoils in cascade. 


ANALYSIS 


The analysis is concerned with an integral equation 
which arises in the problem of an oscillating airfoil be- 
tween walls in two-dimensional subsonic flow. The 
integral equation relates a known distribution of ver- 
tical velocity to an unknown distribution of lift and in- 
volves a complex kernel function. In the first section of 
the analysis the integral equation is presented. In the 
second section the kernel function is examined and re- 
duced to a form suitable for use in calculations. The 
main result of the analysis may be found in Eqs. (8), 
together with the definitions given in Eqs. (9) and (10). 


Discussion of the Integral Equation 


The derivation of the integral equation for an airfoil 
oscillating in the presence of walls has been performed, 


JANUARY, 1955 


for the case of subsonic flow, by Runyan and Watkins. 
In deriving this integral equation one of the boundary 
conditions to be observed is that the vertical induced 
velocity at the walls must be zero. This condition 
can be met by the method of images as has been done 
for the stationary case. The effect of one wall at a 
distance of y = /7/2 on a doublet placed at y = 0 may 
be found by placing an equivalent negative doublet at 
y =H. Itis apparent from the symmetry of the sys 
tem that the vertical velocity at y = ///2 induced by 
the added doublet will be equal and opposite to the ve 
locity induced by the doublet placed at y = 0. Thus 
the boundary condition for the effect of a wall is satis 
fied. To form two walls, by analogy with the steady 
case, it is necessary to consider an infinite row of 
doublets; that is, to add doublets at y = +n// where 
eae | Ds ss ee 

The use of this method of images leads to the follow 
ing expression for the vertical velocity of an oscillating 
airfoil at a point y = 0, midway between two walls 


wh : : 
w(x) = F lim L(x) [A(M, 2) + 
1 


pU~ yO 


K(M, 2, H)]dxo (1 


where w(x) is the known vertical velocity (or known 
motion of the wing) and L(x») is the unknown lift dis- 
tribution or the local strength of a distribution of os 
cillating pressure doublets. The functions within the 
brackets comprise the kernel function of the integral 
equation and appear formally as 


- l ‘ 
aCe, s) = lim e 
DPR y—+0 
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= P+ ro ni) dé 
The first function A(M, z) corresponds to the kernel 
for the free-air condition as given by Possio.® The 
second function K(M, z, HH), containing the infinite 
summation, is the additional part of the kernel arising 
from the effect of the walls. Physically, a kernel func 
tion represents the contribution to the vertical velocity 
at a field point due to a pulsating pressure doublet of 
unit strength located at any other point in the field. 
For the particular case represented by Eq. (2) the 
kernel function gives the vertical velocity in the plane 
of a wing located in the center of the tunnel. The ex- 
pression K(M, z) gives.the downwash of a doublet in 
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the plane of the wing while the expression K(M, z, 7) ing the unknown load, or lift distribution, which pro 
gives the downwash due to the system of images. duces a known vertical velocity. 
With the problem formulated as an integral equation 
(Eq. 1), it may be of interest to note an analogy to the Reduction of the Kernel Function 
structural beam problem wherein the deflection of a = : . ; ; A 
; . [he integrals contained in the expressions for the 
point on the beam can be expressed as the product of a 
load and an influence coefficient giving the deflection 
per unit load. In the present problem the deflection 


kernel function in Eq. (2) are improper since they have 
an infinite limit and also since at certain points the 


; mae i integrands become singular. This section is concerned 
of the point is replaced by the vertical velocity w(x), the oe : - ; : 
. : ; with the reduction of these integrals to a form more 
structural load by the aerodynamic load L(x»), and the ( : 
att : : amenable for computation. 
influence coefficient by the kernel function. In place 


of seeking a deflection produced by a known load, how- Making use of the fact that the Hankel functions in 
ever, the problem here is the more difficult one of find- Eq. (2) satisfy the identity 
O° EP Be pr = 8°O" . [Mk \WI°k Mk 
- — VE + py*)] = — — — 3°) ~ I] Vie + By 3) 
Oy" B- O& pb- 3 re] 


there is obtained for the downwash 
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The integrals of Eq. (4) which contain partial derivatives of Hankel functions can be integrated twice by parts 


to obtain 
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The last integral of Eq. (5) may be written in two parts as 
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The first integral on the right side of equation (6) will be left temporarily in integral form and will be treated 
in a following section (see evaluation of S; following Eq. (13)). 

The second integral on the right side of Eq. (6) has not been integrated in closed form; however, in wind tunnel 
problems, it can be handled conveniently by approximate methods. (An alternative means of treating this in 
tegral, which avoids the approximation but is somewhat more tedious, will be indicated in a following section.) 
A practical assumption which is often made in the analysis of the effect of wind tunnel walls is that the tunnel 
height is considered large compared to the wing chord. With this assumption the argument of the Hankel func- 
tion in Eq. (6) can be written as (in the limit as y > 0) 


| 


Mk _;; = - Mk i +} Mk 
— Vi + B2(nH)? = — pnH +1 nil 
3° B? \6nlT 8 
provided that §/@nH < 1. 
This approximation implies that the airfoil images, and in particular the closest image m = 1, are a sufficient 


distance from the airfoil so that the actual distance V £ + 68°(n//)* may be replaced by the vertical distance 6n// 
of the image above the airfoil. Of course, this approximation does not hold for Mach Numbers close to or equal 
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to unity. 
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With this approximation the second integral of Eq. (4) can be expressed 
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and these expressions may be used to express Eq. (4) as 
wb ' 
w(x) = = L(xo) [A(M, 2) + K(M, 2, H)] dx (8) 
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Eq. (8) together with the definition of Eqs. (9) and (10), 
permits the determination of the effect of tunnel walls 
on a lift distribution L(@)), for a given downwash distri- 
bution w(x). In the next section a method of solution 
of Eq. (8) will be presented. 


METHOD OF SOLUTION 


A method of using Eq. (S) to determine the aero- 
dynamic forces on a wing oscillating in the presence of 
plane walls is now discussed. The method under con 
sideration is one of collocation, similar to that used by 
Possio® and Frazer’ for the case of no walls. The ap- 
proach involves the assumption of an appropriate series 
expression for the lift distribution, substitution of this 
series in the integral equation for the downwash, and 
calculation of the downwash at arbitrarily selected 
points on the chord (control points). Thus Eq. (S) is 
reduced to a set of simultaneous equations the un- 
knowns of which are the coefficients of the assumed ex- 


pression for the loading. 


The Expression for the Loading 

The expression which is assumed for the lift distribu- 
tion is a trigonometric series expansion which satisfies 
the Kutta condition at the trailing edge and which has 
the proper type of infinity at the leading edge. This 


expression is 


V (x — xo)? + B2(nH)? 


L(x) M% — 
. = Ao cot a > > A, sin 70) = L(@) (11) 

pl : Z n=1 
where x» = —cos 6 and the A,,’s are unknown coefti- 


cients to be determined in accordance with the down- 
wash w(x), which is known from the motion of the wing. 
It is desirable to rewrite Eq. (8) in terms of the new 


variable 4, as follows 


re 
w(x) = Uk| } L (09) AK4M, 2) sin 0) dO) + 
Le 0 


{ L(@)K(M, z, I) sin od (12) 
0 


The first integral on the right of Eq. (12) is the in- 
tegral expression first derived by Possio® for the con- 
dition of no walls. Its solution has been treated by 
several investigators and will not be discussed here. It 
can be expressed entirely in terms of the unknown co- 
efficients A, of Eq. (11). In order to evaluate the 
second integral on the right of Eq. (12) it is necessary 
to deal further with the kernel function. 


Further Treatment of the Kernal Function Associated 
with the Walls 


The kernel K(M, z, 1) is the sum of four infinite 
series which can be written as 
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a eee wr oe ; and Chien.’ Introducing the variables p and e, 
) K(M, 2, IT) = (C48; a Coa$o + C383 + C4 (1:3) e 
28 where 
where the .S,,’s denote the indicated infinite summations MeIT x-— x, 
of Eq. (10) and the C,,’s the respective multipliers. a on8 and aes 3H 
(7) Series S; and S, of Eq. (13) may be put in a more 


rapidly convergent form according to Infeld, Smith, — the series 5; and S. can be written as 
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where Euler’s constant y = 0.577215. Series S; may be evaluated by utilizing the expression for S, |Eq. (14 
and integrating the resulting expression to obtain 
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” number of control points and equated to the expression 


With the aid of these series S; to Sy, the kernel A(\/, for the downwash. 
; z, [1) may be evaluated. The second integral on the The expression relating the downwash to the motion 
ite right of Eq. (12) may be expressed in terms of the un- of a wing translating (#) and pitching (a) about an 


known coefficients A, upon substitution of the assumed axis located at x, is 
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wx) =h+ Uat d(x — xa (18) wll 
‘ = mB(2n — 1) (n = 1,2,3,...) (21) 
a 


or, with the assumption of harmonic motion 


+ [1 + itk(x — x,y) Ja (19) 
) 
Eq. (19) is used to calculate w(x) for values of x appro- 
priate to each of the selected control points. <A set of 
simultaneous equations can then be written, the number 
of which corresponds to the number of control points 
employed and (conveniently) to the number of terms 
retained in the series for L(@)). The unknown coeffi- 
cients may now be determined by solving these simulta- 
neous equations. The total lift and moment about 
the midchord are given in terms of the coefficients 


A, through the relations 


—-L 1 -- - 
mobU? 2\°°' 2°" 


M.. l l 
Cao Ay + A» 
apb? S 2 


The Effect of the Number of Control Points Considered 


(20) 


In a series of calculations performed elsewhere, an 
investigation was made of the number of terms of the 
series for the lift distribution [Eq. (11)] and thus the 
number of control points required to obtain satis- 
factory accuracy. Calculations were performed for a 
particular case by increasing the number of control 
points until the solutions were in reasonable agree- 
ment. For the case considered three terms of the 
series for the lift and three control points at the 
1/s-, 1/o-, and */s-chord positions gave satisfactory 
The consideration of two additional control 
trailing edges together 


results. 
points at the leading and 
with two additional terms of the lift series made no 
significant change in the results. For high values of 
the reduced frequency parameter, k, the use of addi- 
tional control points might be necessary. 

The procedure just discussed involves consideration 
of a continuous distribution of pressure doublets over 
the chord. Calculations requiring much less comput- 
ing can be made by considering the chordwise loading 
to be concentrated in a single doublet located at the 
quarter chord and satisfying the downwash at the */, 
chord. In the case of the lift, this approach has been 
found to give fairly good agreement with the results of 
the more elaborate calculations except in the vicinity 
of the resonant frequency. 


DISCUSSION OF THE RESONANCE PHENOMENON 


By examination of Eqs. (14) and (15) it may be seen 
that the series S; and S: become infinite when 


tp? = (2n — 1)? 


or where 


(77 will now be regarded as having dimensions of feet 
for the remainder of the text and in the figures.) At 
these critical values of the frequency parameter the ex 
pression for the kernel A(MV, 2, /7) [Eq. (13)] becomes 
infinite for all values of x. This implies theoretically 
the existence of infinite velocities at the center of the 
These velocities occur in a plane perpendicular 
Physically, this 


tunnel. 
to the direction of the main stream. 
represents a resonance phenomenon and corresponds 
to the well-known case of a simple undamped spring 
mass system where at the resonance frequency the 
equation shows an infinite deflection of the mass. In 
the spring-mass system infinite deflections, of course, 
cannot be obtained; however, the resonant frequency 
is accurately predicted, in spite of the existence of vari 
ous damping factors. Correspondingly, in the wind 
tunnel problem, the infinite values of w(x) would never 
be realized. Under actual conditions factors such as 
finite tunnel length, transmission through walls, turbu 
lence of the flow, and so forth that give rise to damping, 
make pure resonance unobtainable. Even with damp 
ing present, however, resonant frequencies would exist 
at values of w///a in the neighborhood of those given in 
Eq. (21). Experiments have indicated that significant 
resonant effects do, in fact, exist. 

It can be shown also that as the resonant condition is 
approached the loading must approach zero since the 
product of the loading, 1(@)) and the kernel, A (.V, z, /7) 
must approach the known finite value w(x) which rep- 
resents the actual motion of the wing. As A(M, gz, /7/) 
becomes infinite this product can remain finite only if 
L(6)) approaches zero. 

The fundamental or smallest critical values of w///a, 
corresponding to m = | in Eq. (21) are shown plotted 
as functions of Mach Number J/ in Fig. 1. Eq. (21) 
and Fig. 1 show that finite values of the critical fre- 
quency exist for the condition 17 = 0, LU’ = 0, anda # 
©. These conditions correspond to a compressible 
fluid at zero velocity in the tunnel. As the Mach Num 
ber is increased the critical frequency parameter de 
creases rapidly and becomes zero at a Mach Number of 
unity. 

THREE-DIMENSIONAL EFFECTS 


SOME REMARKS ON 


As is well known, the analytical treatment of the 
airfoil of finite span in free space is much more difficult 
than that for the two-dimensional airfoil. Similarly, 
in the consideration of the effect of walls the three- 
dimensional problem is much more involved than the 
corresponding two-dimensional case. The added diffi 
culty arises from the fact that, whereas in two-dimen 
sional flow a single line of images extending to infinity 
is sufficient to satisfy the boundary condition, in the 
three-dimensional problem the images must be dis- 


tributed over the entire plane. It is possible, how 
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ever, to obtain the resonant frequencies for the three- 
dimensional wall problem from simple considerations. 
By use of the method of images, as was done in the 
two-dimensional case, a condition of resonance can be 
shown for the closed rectangular tunnel. For the rec- 
tangular tunnel the fundamental resonant frequency is 


given by the following equation 
wll rae] 
a Vid¢r 


where r is the ratio of the tunnel width to the tunnel 
height. Note that this equation is identical to that 
for the fundamental mode of the two-dimensional 
tunnel [Eq. (21)] except that 7 is replaced by H 
Vi+ 7°, which corresponds to the length of the 
diagonal of the rectangular section. 

A resonance can also be demonstrated for the closed 
circular tunnel. The nature of the boundary value 
problem, for this case, makes it desirable to separate 
variables, so that the governing partial differential 
equation can be reduced to Bessel’s equation. The 
resonant frequencies are then found as the roots of the 


equation 
J ,’(wD/2aB) = 0 ion & 2. fo...) 


or 
wD/2a = p,B 


where D is the diameter of the tunnel and p,, is the root 


of the equation 
J,."(p.) = 0 


Values of p, for the first several modes are p, = 1.84, 
B00, 37,..0%< 

It has been shown that the consideration of three- 
dimensional rectangular or round tunnels leads to 
different values of the critical frequency from those 
obtained in the two-dimensional case but does not do 
away with the possibility of resonance. For example, 
the fundamental frequency for a circular tunnel having 
a diameter equal to the height of a two-dimensional 
tunnel is 1.17 times the frequency of the two-dimen- 
sional tunnel whereas the frequency of a square tunnel 
is 0.707 times that of the two-dimensional tunnel hav- 
ing the same height. By treatments similar to those 
employed herein, it can apparently be shown also that 
under idealized conditions resonance can occur in tun- 
nels of nonuniform cross section (expanding or contract- 
ing section) or even in open tunnels or jets. 


APPLICATION TO SPECIFIC EXAMPLES 


The foregoing theory and calculation procedure has 
been applied to a number of examples. The investi- 
gation has treated three aspects of the wall effects prob- 
lem: (1) calculations have been made to determine the 
lift on an oscillating wing for various Mach Numbers 
at constant wall distance; (2) the effect of varying the 
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distance between the walls is studied; and (3) aerody- 
namic lift and moment coefficients obtained for an air- 
foil oscillating between walls have been employed in an 
idealized flutter analysis and results have been com- 
pared to those for an airfoil in a free stream. 


Effect of Variation in Mach Number at Constant Wall 

Distance 

Application has been made to the particular case of a 
wing oscillating in pitch about its midchord. Calcula- 
tions have been performed for the height-chord ratio 
A = 3.80, and for tunnel speeds corresponding to 17 = 
0.3, 0.5, 0.7, and 0.8. Results of these calculations are 
shown in Fig. 2. The vertical axis is the ratio of the 
magnitude of the lift | 1} on a wing between walls to the 
lift | L’ 
indicates no wall effects. 
plotted against the frequency parameter w///a for the 
four values of Mach Number mentioned previously. 
The curve for each individual Mach Number covers a 


on a wing in free air so that a value of unity 
This ratio of the lifts is 


range of values of w///a from zero to the value for the 
fundamental resonant frequency. The critical value 
corresponds to the value shown on the previous figure 
and is again seen to decrease as the Mach Number in- 
creases. 

Note that at a particular value of the frequency 
parameter a progressively larger effect of the walls is 
indicated as the Mach Number increases. At all 
Mach Numbers the lift is reduced to zero at the res 
onant frequency. The dip in the curves against fre- 
quency ratio, which seems to be characteristic of the 
low Mach Number cases, gradually disappears as the 
Mach Number is increased. Initial calculations for the 
moment on the wing indicate that curves of essentially 
the same form will be obtained. 

For flutter and for other problems of nonstationary 
flow, knowledge is required not only of the magnitudes 
of the forces involved but also of the phase relationship 
of the forces to the motion of the wing. A phase rela- 
tionship between lift force and position of the wing is 
shown in Fig. 3 for the cases just discussed. The 
phase relationship is shown as the difference (¢ — $’) 
between the phase angle for a wing in a tunnel and that 
for a wing in free air and is plotted against the fre- 
quency parameter, w///a, for four values of the Mach 
The effect of the tunnel walls appears as a 
As in the 


Number. 
deviation from zero of the value of ¢ — @’. 
case of the magnitude of the lift, the effect of the walls 
on the phase angle increases as the resonant frequency 
of the tunnel is approached and also as the Mach 
Number is increased. Although it is not shown by this 
form of presentation, it should be stated that the phase 
angle, @, for the wing in the tunnel attains a value of 


90° at the resonant frequency. 


Effect of Varying the Distance Between the Walls 


The effect on the lift force ratio of changing the dis- 
tance between the walls is illustrated in Fig. 4. The 
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AIR FORCES ON A WING 
results presented in the two previous figures have been 
based on the consideration of a distribution of pres- 
sure doublets over the chord of the airfoil and on satis- 
fying the downwash condition at three chordwise sta- 
tions. Results for Fig. 4 have been obtained by the 
simplified procedure of concentrating the loading in a 
single doublet at the quarter chord and satisfying the 
downwash condition at the three-quarter chord. As 
stated previously this approach gives fairly good agree- 
ment with the results of the more elaborate procedure 
except near the resonant frequency. 


Calculations have again been made for an airfoil 
oscillating in pitch about its mid-chord for values of the 
) 


tunnel height-chord ratio \ of 3.80, 8, and 16, and for 
M 0.3 and 0.8. 


is plotted as ordinate and the 
Plots for 


The lift ratio | L/L’ 
reduced frequency k = bw/l’ as abscissa. 
both Mach 
for ease of comparison. 
the effect of reducing the Mach Number is to reduce the 


Numbers are made to the same _ scale 
It is again apparent that 


effect of the tunnel walls and to raise the value of the 
critical frequency at which resonance can occur for a 
given tunnel. For example, for \ = 3.80, at 17 = 0.8 
the critical k is 0.31, whereas for J = 0.8, the critical 
k is increased to 1.51. Also, as was to be expected, in- 
creasing the height of the tunnel has a marked effect 
in reducing the influence of the tunnel walls for most of 
the k range. However, the critical frequency is re- 
duced for increased tunnel height so that in the large 
tunnel the range of k below the fundamental resonance 
disad- 


becomes smaller. This would seem to be a 


vantage of the larger tunnels. However, if the second 
branch of the curve for \ = 16 at 1/ = 0.3 is examined, 
it is seen that, for frequencies between the first and 
second resonance points, the effect of the walls on the 
magnitude of the lift is no greater than for frequencies 
below first resonance. Note also that the approach 
to resonance is quite abrupt, so that only a small range 
of frequencies very close to resonance would be critical, 
thus permitting tests to be conducted between the criti- 


cal frequencies. 


As a matter of interest some results of Reissner* 
(lJ = 0, A= 3.14) are compared to results of the pres- 
ent analysis (AJ = 0.3, = 3.80) on Fig. 5. The curve 
for lJ = O significantly duplicates the rather large wall 
effect at low values of k which is noted for the JJ = 0.5 
result. [At low values of k the curves for the two 
different Mach Numbers are coincident probably be- 
cause of the fact that the slightly lower height-chord 
ratio (A = 3.14) used by Reissner counteracts the effect 
of decrease in Mach Number.| For values of k greater 
than 0.5 the curves for the two Mach Numbers separate 
and the lift force ratio at .1J = O approaches unity, ex- 
hibiting no effects of resonance, since the resonance 


phenomenon arises only from the effects of compressi- 
bility. 


CSCILLATING 


WALLS 1%) 


BETWEEN TWO 


Flutter of an Airfoil Between Walls or in Cascade 


In addition to the effect of the walls on the forces on 
an airfoil, it may also be of interest to examine the 
effect of the walls on the calculated flutter speed of an 
airfoil. The results of such calculations can also be 
interpreted in a limited sense as the flutter of airfoils 
in cascade, since the image system used to satisfy the 
boundary conditions for tunnel walls actually repre 
sent airfoils oscillating 1SO° out of phase with each 
other. Lilley'’ has found experimentally that it is 
possible for some types of flutter in cascade to occur in 
this manner. Of course, these results are based on the 
concept of potential flow at small angle of attack but 
may provide a framework with which to study the high 
angle of attack stall type of flutter which is frequently 
encountered in compressor blades. 

For these calculations, use has been made of the air 
0.7, and \ = 3.80 and the 


The ordinate is the flutter 


forces already shown for .\/ 
results are shown in Fig. 6. 
speed coefficient l'/bw, and the abscissa is the fre- 
quency ratio w),/w,. 

Curves are presented for a wing in a tunnel (or 
cascade) and in free air for two different values of the 
*nass parameter |/x. The wing characteristics are shown 
on the figure and apply in this case to a two-dimensional 
wing mounted on springs having a frequency in rotation 
of w, and a frequency in bending of w,. The curve for 
1/k = 500 corresponds to conditions in a compressor, 
whereas the curve for |/k 50 corresponds to values 
of the mass parameter more normal for aircraft wings. 
For 1/« = 500, the effect of the presence of airfoils 
in the cascade seems relatively small for the case con 
However, the gap-chord ratio of 4/1 is large 
Further cal 


sidered. 
compared to the more normal ratio of 1/1. 
culations are needed in order to more completely study 
For the condition representing a wing 
50, the effect of the walls is 


the problem.* 
in a tunnel, 1.e., for 1/«= 
larger than at the higher value of 1/k. 
is of the order of 


A practical 
value of the frequency ratio w,/w 
0.4 to 0.5. In this range the effect of the tunnel wall 
changes the flutter speed 12-15 per cent. Note that 
for free air, the curve is continuous for the complete 


a 


range, whereas the curves for a wing in a tunnel 
This break in the curve is due to 
In this case the 
.< 1.1. The 
two branches of the tunnel wall curve are associated 
with flutter frequencies below the fundamental reso 


Wy / Wa 
are not complete. 
presence of the resonant condition. 


resonance occurs in the range 0.7 < w/w 


nant frequency of the tunnel. 


* For application to actual cascade configurations, or in wind 
tunnel problems where the gap to chord ratio is small, the ap 
proximation employed in integrating Eq. (7) becomes less valid. 
It is possible to avoid the use of the approximation by writing 
the integral of Eq. (7) in a form which is identical to that of Eqs 
(16a) and (16b) with the exception of the upper limit. The in 
tegral containing the Hankel function can be evaluated by em 
ploying the tables of Schwarz. The second integral, containing 
only an exponential term, can be integrated in closed form, as was 
done to obtain Eq. (16e 





50 JOURNAL OF THE 


CONCLUDING REMARKS 


This paper deals with the forces on an oscillating air- 
foil between plane walls in subsonic compressible flow. 
Of primary concern has been the presentation of a 
method for calculating these forces and the application 
of the method to various examples. 

A resonant phenomenon, involving a transverse os- 
cillation of the moving air between the walls is dis- 
cussed. The conditions under which resonance can 
occur, previously pointed out by Runyan and Watkins 
for the case of two-dimensional walls, are given for 
three-dimensional tunnels of rectangular and circular 
cross section. It is shown that at resonance the load 
on the oscillating wing approaches zero, and conse- 
quently experimental testing in tunnels should be done 
under conditions as far removed from resonance as 
possible. 

Numerical calculations have been made to determine 
the forces on a wing at various Mach numbers and for 
various wall heights. Frequencies of oscillation from 
zero to the fundamental resonant frequency have been 
considered. As was to be expected, the effect of the 
walls on the lift ratio and the phase angles was shown 
to be more pronounced as the Mach Number was in- 
creased and, at high Mach Numbers was found to be 
large even at frequencies well removed from resonance. 
The effect of increasing the height of the tunnel was to 
greatly diminish the effect of the walls but also to lower 
the resonant frequency. 

A few calculations were made at a Mach Number 
of 0.7 to determine the effect of the tunnel walls on the 
flutter speed of a simple configuration and to show the 
effect on flutter of airfoils in cascade. Some represen- 
tative calculations for a wing in a tunnel showed a re- 
duction in flutter speed of about 10 to 15 per cent from 
the values calculated for free air. For the flutter of 
airfoils in cascade, at a value of 1/x representative 
of compressor blades, it was found that the effect of the 
presence of an infinite array of airfoils oscillating 180° 
out of phase was not very great. However, the gap to 
chord ratio used was not realistic for actual configu- 
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rations so that further work is needed before any con- 
clusions can be made. 

The calculated results are naturally affected by such 
assumptions of the theory as the absence of damping 
(which would arise from finite tunnel length, trans- 
mission through walls, turbulence of the flow, and so 
forth) so that in actual tests the wall effects would be 
somewhat less pronounced, but could, under certain 
conditions, still be large. The investigation just dis- 
cussed permits the regions where wall effects may be 
greatest to be clearly defined and perhaps provides a 
means for evaluating experimentally measured data. 
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Secondary Flows in Cascades 
of ‘Twisted Blades 


FREDRIC F. EHRICH* 


Westinghouse Electric Corporation 


SUMMARY 


§ Squire and Winter's? definitive paper on secondary flow treats 
the situation of a nonuniform velocity profile through a cascade 
of blades whose properties—i.e., geometry—are uniform with 
The analysis focuses on the flow in a single passageway 


has treated the special case of uniform 


height. 
between blades. Tsien 
a cascade of blades whose properties vary as a 


He uses a classical Prandtl lifting line 


flow through 
function of blade height 
inalysis 

Using a conceptual model similar to that of Squire and Winter, 
. general analysis is carried out to determine the perturbation 
or secondary flow resulting from a nonuniform velocity profile 
passing through a cascade of blades whose properties vary as a 
function of blade height. The general result reduces exactly to 
the special form of Squire and Winter and also gives a result 
quantitatively and qualitatively equivalent to that of Tsien 

Experimental measurements in a cascade give results indicating 
the pertinence of the analysis except in the regions where viscous 


stresses invalidate the assumptions 


SYMBOLS 


= Eigen value 
1 = coefficient of the homogenecus solution 
C a constant 
= chord length 
coefficient for comparison of passage and lifting line 
analyses 
a particular solution of the partial differential equation 
R = slope of the (lift )-(angle of attack) curve 
l = blade length 
L = lift 
} = harmonic coefficient 
static pressure 
= total pressure 
= pitch 
l = velocity in the direction of flow 
= velocity component in the horizontal direction 
= velocity component in the vertical direction 
\ horizontal coordinate 
vertical coordinate 
Y = function of x 
a function of y 
= coordinate in the direction of the main flow 
a = velocity angle component in the vertical direction 
blade circulation 
= harmonic component of blade circulation 
j = deviation of inlet angle from the mean 
6 = deviation of exit angle from the mean 
turning angle 
harmonic component of turning angle 
6 = angular direction of a streamline 
= vorticity component along a streamline 


Presented at the Aerodynamics Session, Twenty-Second An- 
nual Meeting, IAS, New York, January 25-29, 1954 
* Mechanical Engineer, Advanced Development Section, Avi- 


ition Gas Turbine Division 


p = fluid density 

rr) = angle between the normal toa Bernoulli surface and the 
principal normal of a streamline 

Vv, = stream function of the flow in the exit plane 

Vv, = stream function of the two-dimensional or design flow 

y = stream function of the perturbation or secondary flow 


See also Fig ] 


Subscripts 


P = particular solution 
H = homogeneous solution 
0 = reference 


INTRODUCTION 


has derived an expression for the 
the 
streamline executes an arbitrary path 


a 
change of 
streamline as the 
with no dissipation in an incompressible medium: 


vorticity component along a 


9 


V(P, p)| sin ¢@ da/Ul? (1 
/ 1 

An effective procedure for an analytic description of 
the secondary flow in cascades is to make the following 
simplifying assumptions: 

(1) The flow is assumed to be incompressible and 
inviscid. 

(2) It is assumed that there is a velocity variation in 
This 


assumption is accurate for axial flow turbomachinery 


the height direction only entering the cascade. 


since alternating rows are rotating with respect to one 
another and wakes of the preceding blade row are seen 
as velocity fluctuations in time rather than in space. 

(3) The radius of curvature of the blade passage is 
assumed large enough so that the pressure gradients 
across the blade passage are negligible. This seems rea 
sonable in axial compressor work but is questionable 
in turbines where very large turnings are encountered. 

(4) Pressure variations in the height direction are 
neglected. This limits our attention to cascades of par 
allel blades or annular cascades with hub-tip ratios very 
close to unity. The blade twists are considered in their 
relation to secondary flows and not to centrifugal equilib 
riums. 

(5) The angle of turning is assumed small enough 
so that the plane Bernoulli surfaces are not appreciably 
distorted. This is probably the grossest of the as- 
sumptions in that it limits the solution to a mere per 
turbation on the main flow. Hence, the solution is 
more indicative of the trends existent rather than the 


actual resultant conditions. , 
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SECONDARY FLOWS IN CASCADES OF TWISTED BLADES 53 


ANALYSIS 


Considering a plane of the blade passage as in Fig. 


2, the change in vorticity is found to be: 
~ § = —2(0U/dy)e (2 


If the blade passage is twisted, the inlet and exit 
angles at an arbitrary height point will differ slightly 
from the angles at the mean plane. For simplicity, 
let us consider the inlet and exit section at two stations 
where the blade passage section is not appreciably dis 
torted from a rectangle, as in Fig. 3. 

Assuming the existence of plane flows at the inlet, 
examination of the section ‘‘A—A”’ at inlet reveals the 
existence of a set of cross flows as in Fig. 3. 

These flows have the velocity magnitude: 

u(x, y) = Ub, (3a) 
v(x, vy) = O (Sb) 


so that the vorticity of the flow at the inlet will have a 


magnitude: 


1 = Ov/Ox — Ou/Oy = —O(U5,)/Oy (4 


Str 


Hence the net vorticity of the flow at exit will be: 


= —2«(0U/dy) — O(U6,)/Oy (& 


&2 
From the equation of continuity for the section 
“B-B” at exit, we may define a stream function for the 
flow in this plane neglecting changes in the main direc- 


tion of flow such that: 
u = OW,/Oy (6a) 
v= —OWV,/Ox (6b 
from which: 


= —Ov/Ox + Ou/Oy = VY, = 
2€(OU/Oy) + O(U6;)/Ov (7 


The boundary conditions of this equation are deter- 
mined by the flows normal to the boundaries. At the 
hub and tip walls, the flow normal to the boundaries 
is zero, so that: 


v(x, 0) = 0 W.(x,0) = C (Sa) 


v(x, 1) = 0 V(x, 2D =C (Sb) 


Along the sides of the passage, we will again find 
cross flows due to the twist of the passage as in Fig. 3. 


u(+s/2, y) = Ub. W,(+5/2, y) = { Ube dy (8c) 


Eq. (7) and boundary conditions of Eq. (8) com- 
pletely determine the flows in the exit section ‘“B-—B”’; 
but in order to interpret the secondary flow more 
clearly, it is convenient to consider it as the summa- 
tion of two flows—first, the two-dimensional or ‘‘de- 
sign’’ cross flow and secondly, the additional disturb- 
ance or ‘‘secondary”’ flow on this two-dimensional flow. 


Specifically, the desired two-dimensional design flow 


at the exit will be: 


v(x, ) 9b 


as in Fig. 4. 
This flow can also be described by the equation and 


boundary conditions: 


vw: O( Lb.) Ov LO 
v(x, QO) = @ W(x, O ( lla 
\ 0 Wo! a llb 
u( +s/2. 4 l6 Wo(+5/2, 4 


Defining our perturbation stream function as 


y VY, — 12 
which is the difference between the actual and two 
dimensional flow in section ‘‘B-B,’’ we may subtract 
the set of Eqs. (10) and (11) from Eqs. (7) and (8) giv 
ing: 


V°y = 2€(O0l’ Oy) + O(U6;) OV — O(L/b2)/O Ls 


vx, UG) = G V(x, 0 0 (l4a 
v(x, /) = 0 (x, / 0 l4b 
u(+s/2, y) = 0 ¥y(+s/2, y) = 0 (l4e 


The angle of turning at an arbitrary height point as 
shown in Fig. 3, is: 


€ = «o + 6, — bo lo 
So that Eq. (13) may be rewritten : 


V°y = 2e(0U’/Oy) + U(Oe/Ov) + (6; — 62) X 
(OU’/oy) (16 


Since by assumption the third term is small com 


pared to the first, we get: 
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or finally: 
s 1 O(eU?) 
7 l oy — 
v(x, 0) = 0 (19a) 
y(x,/) = 0 (19b) 
¥(+s/2, y) = 0 (19¢c) 


ANALYTIC COMPARISONS 


In the case of no variation in the turning angle in 
the height direction we find, of course, that the solution 
of Eq. (18) is identical to the results of Squire and 
Winter? giving: 


€ = & (20a) 


V*y = 26 (OU /Oy) (20b) 


In the case of a uniform inlet velocity profile, the 
secondary flow is described by: 

U = Uy (21a) 

V*y = Up(dc/Oy) (21b) 

Solution for the secondary velocity distribution is 

accomplished by first solving Eq. (21) for the stream 

function by means of separation of variables (Appendix 

I). This gives: 
y= Lj1- 

1 nw 


n cosh (n7rs/2/) 


cosh (nx //) | Ul 


4 nTy 
€, Sin ] (22) 
F >] 
2 nTy : 
é, = € COS dy (23) 
| a are l r 


From the definition of the stream function, the vertical 


where: 


component of the velocity ts: 


oe Oy U, 3 (= (mmx | 
1 


Ox cosh (7s /2/) 


, nTry 4 
€, Sin (24) 
l 


If we consider the section under study to be imme- 
diately before the trailing edge, there will be a discon- 
tinuity across the blade in the velocity component in 
the height direction. In the section immediately after 
the trailing edges, this discontinuity will, of course, 
persist and hence be exactly equivalent to a vortex 
sheet. The strength of this vortex sheet may be calcu- 
lated by an examination of the velocity discontinuity 
as in Fig. 5 for an arbitrary secondary flow. Note that 
we may study the perturbation stream-function alone 
since the two-dimensional flow which we have sub- 
tracted out has no velocity components in the y direc- 
tion. 

The circulation for a short length of this vortex sheet 


is then: 
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dV = v2 dy — v, dy (25 


Due to the symmetric nature of the secondary flow 
analysis, the height-wise velocity component on one 
side of the trailing edge is equal in magnitude but op- 
posite in sign to the velocity component on the other 
side. Thus, the strength of the vortex sheet or circu- 
lation per unit length is equal to: 


(26a 


dV /dy = —2v],~,/2 


so that, from Eq. (24): 


—— nars\ . nTry ' 
—2U 2. «, tanh sin (26b 
n=1 2] l 


Since we are dealing with the special case of uniform 
inlet velocity or potential flow, the problem is also amen- 
able to treatment by the classical Prandtl lifting line 
This introduces the concept of shed circu- 


dv dy = 


analysis. 
lation to calculate the magnitude of the velocity dis- 
continuity across the trailing edge. 

Tsien® has treated the case of a uniform flow through 
a cascade of twisted blades. Using his results, one may 
calculate the strength of the vortex sheet shed at the 


trailing edge (Appendix II). This gives: 


ne NTs ny 
dV/dy = —2Uy > €, tanh ( = sin ( *) Goa 


n=1 l 


(27) 


where: 
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SECONDARY FLOWS 


l 
l tanh (7s /2/) 


oo (ns /2/) 


It is seen that, except for the coefficient C,, this 
completely different approach to the problem gives an 
answer identical to the result of Eq. (26b) from the 
The coefficients C,, plotted in Fig. 
to (2.000) but through- 


passage analysis. 
6 varies in value from (0.667) 
out much of the important part of its range its value 
is roughly unity so that both analyses give similar 
quantitative results. A sample calculation of a set 
of cascade blades whose turning varies linearly with 


span is given in Fig. 7. 


DESCRIPTION OF EXPERIMENTAL WORK 


Considering the governing equation of secondary 


flows in cascades of twisted blades: 
V*y = 2€(0U’/Oy) + U(dce/dy) (29) 


the first term in the right-hand side of the equation 
represents the secondary flow introduced by develop- 
ment of vorticity in the direction of flow as the fluid is 
turned in the cascade. The second term represents the 
secondary flows introduced in the flow by variation of 
turning angle. The presence of the second term is 
primarily a function of design and offers the oppor- 
tunity of controlling or eliminating the effects of the 


first term as desired. Complete elimination implies: 


1 O(el) 
: = 0 (30a 
l Ov 
or: 
el’? = constant (30b) 


An experimental program was undertaken to determine 
the accuracy in the analytic description of the phe- 
nomena. If it is desired to reduce secondary flows, Eq. 
(30) implies that as LU’ decreases, € must be increased. 
The experimental program involved three steps build- 
ing up to the most general case. 

Case I—A set of blades with nonuniform turning 
was tested in an inlet flow with uniform velocity profile 
to isolate and study the term: [U(0e/Oy)]. 

Case II—A set of blades with uniform turning was 
tested in an inlet flow with a nonuniform velocity pro- 
file to isolate and study the term: [2«(0U//Oy)]. 

Case III—A set of blades with nonuniform turning 
was tested in an inlet flow with a nonuniform velocity 
profile to study the superposition of the two contribu- 


tions. An attempt was made to hold el? constant. 


The Cascade 

Tests were run in a small 3- by 7-in. section cascade 
tunnel at a velocity pressure of 10 in. of water and dis- 
charging into the atmosphere. 
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The blades used were made of rolled stock with the 
following specifications: 
blade section = sharp-nosed compressor section 
blade thickness ratio = 0.11 
blade length (overall) = 3 in. 
blade camber angle = 19.5 
blade chord = 1.441 in. 
pitch = 11/2 in. 
number of blades = 5 
stagger angle (based on tangent chord line of un 


twisted section) = 10 
inlet angle (measured from axial direction) = 14 
turning angle of untwisted section = §S 


Variation in turning angle was obtained by bending 
the rolled stock on a die to give a maximum twist of the 


tangent line with the nose as fulcrum of 10°. 


Velocity Profile 

Since no boundary-layer bleed was available on the 
tunnel, the uniform velocity necessary in Case I was 
approximated by peeling off the boundary layers at the 
inlet to the cascade with thin retaining walls which 
separated them from the main flow. 

A velocity gradient in the inlet profile for Cases I] 
and III was obtained by placing a coarse screen half- 
way into the flow passage at 10 chord lengths upstream 


of the cascade. 


Instrumentation 

Instrumentation was similar to that used by Van Le.‘ 
Velocity angles were measured by means of a cobra 
probe, a pair of 1/32-in. tubes whose ends were beveled 


toa 60° angle. 
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Fic. 7. Comparison of the two analyses for a typical calculation 
of the shed vortex sheet due to nonuniform turning angle 


The total pressure was measured by means of a Kiel 
probe which is relatively insensitive to misalignments 


in orientation. 


Computations 
The stream function is defined by: 
v = —Oy/Oor (31) 


so that the stream function may be calculated from the 
measured velocity component by means of the relation- 


ex 
—_ ey ° (29 
nee, - = | vi dx > 4 
/70 


The velocity component was evaluated by means of the 


ship: 


. 


approximation : 


V = Usina = Ua = Ur — P, p| ce 


The total pressure and vertical velocity angle com- 
ponent were measured directly and the static pressure 
was taken to be approximately equal to atmospheric 
pressure. The quantities were measured over an entire 
passage area on a '/s-in. grid in the plane of the trailing 
edges. 

The integration was approximated from the point-by- 
point data by a simple finite difference summation. 


Thus: 
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v1, const. ~ — { [2(P, — P,) p| adx (3 
J0 


n 


Wy=const. * — DY [2(P: — Ps)/p]“an(Ax) (35 


n=O 


or making the quantities nondimensional : 


_—- (=) > (“*) Pin — =| 
Uoeol LJ n=0\ & Po — P, 


In actual practice, the reference level of a had to be 


adjusted as in reference 2 in order to satisfy the zero 
stream function condition on both boundaries. This 
adjustment was of the order of +0.5° which was also 
the expected accuracy of the angle measurement. 


RESULTS 


Taking the inlet velocity profile as measured, and 
the turning angle variation as measured, the right 
hand side of Eq. (29) may be computed for the three 
experimental cases studied An electrical analogue for 
Poisson's equation has been designed and constructed.® 
Using this analogue, the solution to Eq. (29) with these 
three Poisson terms were readily found to +2 per 
cent. This set of analytic solutions may be compared 
to the measured values of the stream function. All of 
these results are shown in Fig. 8. 


DISCUSSION AND CONCLUSIONS 


The comparisons of the experimental and analytic 
results in Fig. 8 show reasonably good agreement in the 
midportion of the passage. A particularly strong 
deviation from the analytic evaluation is in evidence 
at the hub and tip regions in all three cases. This is 
most probably due to the growth of a wall boundary 
layer within the cascade after the original wall boundary 
layer had been separated off at the entrance. The 
directions of rotation of these extra disturbances are 
in accord with this physical concept. This strong 
re-emergence of a secondary flow induced by a wall 
boundary layer would seem to minimize the possible 
gain of boundary-layer bleed in a turbomachine since 
bleeding would probably be necessary at every blade 


row. 
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APPENDIX I) ANALYTIC SOLUTION TO POISSON 


EQUATION 
\ particular integral of the differential equation: 


y(x, 0) = 0 
vix, 1) = 0 Al 


¥(+s/2, y = () 


is seen immediately by inspection to be 


sf / 
il | | ‘o. ypr(x, 0) = 0 A2 


l Oc Oy) dy ay 


v-y = (0c Ov 


vi F(y) = 
r ¥p( +s/2, y) = Fly (A3 
l € =: ie ay . . 
/7 
This particular solution leaves the homogeneous equa 
tion 
vi(x, 0) = 0 
Viva 0 Wal Xx, / = (A4 
Va(+s/2,v) = —Fly 


Taking the solution of Eq. (A4) to be separable® and of 


the form 
Yo = X(x)V(y (AS 


we find by substitution in the partial differential equa 
tion the two ordinary differential equations: 


XxX" —-a@xX =0 AG 
”+a?Vy=0 AZT 
These have solutions: 
Y A, cosh (ax) + Ao» sinh (ax AS 
}° = A; cos (ay) + Ag sin (ay (A9 


Considerations of symmetry imply: 
Ae = G A10 


while the boundary conditions give: 


a NW / nu = LZ, :. oe Al2 


This reduces the homogeneous solution to: 


* nnrx\ , ny 
Wu > pH A,, cosh ( ; ) sin ( ) A13 


On the boundaries x = +s/2: 


= NWS nNwrvV 
Vu F(y) = = A cosh ("7") sin ( ) Al4 


~dF(y) 
= —li(e — &) = 
Ov 
os nw nNqTSs nwy 
z: A, ( 7) cosh ( ) cos ( 2) (Al5 
1=1 ij 2/ l 
Taking 
—. NnqTV 
e= > «, cos v\ | € A116 
n=1 1 / 
we find that: 
Ue 
{ = A173 
(nw /l cosh (7s 2] 


giving 


_] > € cosh Nw. 


nTry 
sin ( ’) AIS) 


Combining the particular and homogeneous equations 


1 nw cosh (mrs 
then gives 


1 WW 


cosh (nx / ny 
~ sin : AI9) 
cosh (mms 2 
or finally: 


‘osh (max / 
y= GD 7 YT x 
oe. 9) 


cosh (ms 21) J 


APPENDIX II -TSIEN’S CALCULATION OF SHED 


CIRCULATION 


Tsien® has studied the situation of uniform flow 
through a cascade of twisted blades by the classical 
lifting line concept. He finds for the circulation dis 
tribution along the blade: 


? Ustck X 
ij T a NTS nry 
€— :  » n y, coth cos 
1] 1 2] l 
= ny 
pl z 7% cos ( 2) \?1) 
n=0 l 
where: 
< {ny 
T= ia ™ cos ( aed AD?) 
0 l 
Taking: 


we may solve for y, and find: 
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€, 
Yn is \ A24 
2 Tn nS 
~ coth 
Lock WU, \ 2/ 
giving for the circulation: 
a é, nTry 
I z cos : A?5 
n=0 2 7 NTS l 
t — coth 
lock Uo 2] 
and for the shed circulation: 
dV dy 205. 2a MX 
Nn l 
NTS 
e, tanh 
2/ . [nr 
sin (A26 
tanh (ms /2/ l 


1/ 


») ? / 
MTC k i 


A more direct equivalence between the lifting line 
analysis and the passage analysis can be had if we notice 
that in the former case the two-dimensional lift of the 


airfoil section is given by: 
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= (pl 2)cR «€ (A27 
while in the latter case: 
L. = puT pU7 es (A28 
so that: 
ck — 2s (A29 


We find then for the final evaluation of the shed cir- 


culation: 


dV dy - 2U)y bf 


nas\ . [ny 
e, tanh ( ) sin ( ”) iC]. (Ago 
2] l 


where 


tanh (n7s/2/ 
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The Effects of Changes in Aircraft Weight on 
Flight Equilibrium 


oan Birman 
ngineer, Systems Engineering Department, The W. L. Maxson 


Corporation, New York 
June 10, 1954 


a 


I 4 RECENT PAPER, Edward S. Rutowski derived a formula 
reference 1, page 195) relating the change in gross weight 
of a jet-powered aircraft to the change in altitude during the 
We have had occasion to make a 
i performance problem, 


cruise portion of the flight 
similar calculation, in connection with 


it a somewhat different formula. In particu 


and have arrived 
lar, Mr. Rutowski neglects the thrust-drag equilibrium require 
ment and assumes that the aircraft maintains constant speed 
whereas our calculations indicate that these approximations (par 
ticularly the former) are not justified 

The problem is as follows: If an aircraft with gross weight 
W is initially in level, unaccelerated flight at a specified altitude 
and speed, what changes will occur in its altitude and speed if 
the weight changes to (HW + AW)? 

The conditions for steady-state equilibrium during flight are 


C,(.S/2)p V’? = W l 
Ca(S/2)9 V2 = : 2 
Cm =U 3 
where 

( = lift coefficient 

C = drag coefficient 

Cm = pitching moment coefficient 

Wo = weight 

7 = thrust 

p = air density 

V = airspeed 

» = wing area 


If there is a decrease in the weight, W, then the lift must also 
decrease, in that restored 
This can be accomplished in either of two ways: 

If it is desired to maintain the initial altitude and speed, 
then the lift coefficient, C,, must be adjusted to rebalance Eq 


order steady-state equilibrium be 


(] However, any changes in angle of attack would unbalance 
Eq. (3), and therefore the pilot must the aircraft by ad 
justing the position of the elevator control so that Eq (3) is 
At the same time, the pilot 


“trim” 


satisfied at the new angle of attack. 
must adjust the throttle position so that the thrust is equal to 
the drag at the new angle of attack 

b) If the pilot does not change the position of the elevator and 
throttle settings, then the altitude and speed must change until 
equilibrium conditions are restored at the initial angle of attack. 
lhe required changes in air density and speed in case (b) can be 
determined from the condition that Eqs. (1) and (2) must be satis 
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fied at the new equilibrium position 


C,(.S/2) (p + Ap) (V + AV)? W AW (4 
JS oe ol 2 P 
Ca (p + Ap | + AV)? = 7 tT Ap + Al oO 
4 op ol 
Combining Eq $4) and (5) with Eqs 1) and (2), neglecting 
second-order terms, we obtain 
(W/p\Ap + (2W/WAV = AW 6 
{((T/p) — (OT /dp)| Ap + [(27/V) — (OT/OV)|AV 0 7 
Simultaneous solution of Eqs. (6) and (7) for Ap and AV gives 
Ae i oz 
Ap VW W ov 
- = - 8 
All 20/7 1 o7 
V Op pov 
1 oT 17 
AV W Op oll 
= + J 
AW 2 o7 1o7 
V Of poVv 


rhe required change in altitude can be found from Eq. (8) an 
the additional! relationship: 


Ah/AW = (Ap/AW) X (Ah/ Ap) 10 


where Ah/Ap is determined from the “standard atmosphere’’ 


relationships 


The quantities Ah/AW and AV/AW were computed from Eq 


8), (9), and (10) for typical combinations of V and h. The 


weight, W, was taken to be 14,000 Ibs., and the wing area, 5S, 
was taken to be 250 sq. ft. The thrust, 7, which is required 
was determined from Eqs. (1 
ind (2), using The 
O07 /Op and OT /OV were determined from the engine performance 
data for the Pratt & Whitney J-42-P-4 engine.? For that engine, 


there is a fixed relationship between the throttle position and the 


it different values of p and V, 


tvpical aerodynamic data. derivative 


engine speed, so that constant throttle setting implies constant 
engine speed. The derivatives 07/0p and O7/OV were there 
fore computed from the curves of available engine thrust plotted 
against speed and altitude, at fixed engine r.p.m. The results of 


the calculations are summarized in Table 1. 


TABLE | 


Knots Ah/AW AV/AV 
Ft. per Lb.) Knots per Lb.) 


Altitude h (Ft.) Speed V 


15,000 340 3.46 6.99 x 10 
480 3.88 12.7 X10 
30,000 340 3.33 8.69 x 1073 
480 3.69 10.4 x 10 
35,000 340 2. 67 6.08 * 1073 
480 3.10 12.8 xX 1073 
40,000 340 2.27 6 18 kK 1078 
480 3.01 i7.6 xX iors 








62 JOURNAL OF THE 


For a basic aircraft weight of 14,000 Ibs., Mr. Rutowski’s 
derivation gives the numerical results 
Ah/AW = —(20,940/W) = —(20,940/14,000) = —1.5ft. per lb. 
AV/AW =0 
for altitudes above 35,000 ft. and all speeds. His formula was 
derived on the assumption that Eq. (1) can be rebalanced after 
a change in weight by an appropriate change in p. He neglected 
the fact that any changes in p will act to unbalance Eq. (2) and 
that both Eqs. (1) and (2) must be satisfied for equilibrium dur 


ing flight. 
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ance Problem, Journal of the Aeronautical Sciences, Vol. 21, No. 3, pp. 187 


195, March, 1954 
2 Specification No. N-1603-C, Pratt & Whitney Aircraft Division, United 


Aircraft Corporation, East Hartford, Conn 


Author’s Reply 


Edward S. Rutowski 
Douglas Aircraft Company, Inc. 
October 29, 1954 


M* BIRMAN claims that the subject paper! neglects the 
thrust-drag equilibrium requirement in deriving a formula 
relating the change in gross weight of a jet-powered aircraft to 
the change in altitude during the cruise portion of the flight 
Her claim, however, is not justified in that the thrust-drag 
equilibrium requirement is implicitly satisfied for a turbojet 
powered aircraft flying at a constant attitude and at a constant 
Mach Number if the thrust is assumed to be directly proportional 
to the atmospheric density. This frequently made assumption, 
which was implicit in the derivation in the subject paper, is rea- 
sonably good for a turbojet at altitudes above 35,000 ft. in the 
isothermal layer since the turbojet is essentially a temperature- 
limited device. The engine characteristics used by Mrs. Birman 
apparently do not fit this assumption well, and hence her results 
differ from those in the subject paper. 

If Mrs. Birman had chosen an engine with the characteristic 
of having a thrust directly proportional to the atmospheric den 
sity at a constant Mach Number above 35,000 ft., her results 
and those in the subject paper would have coincided. This may 
be seen by evaluating Mrs. Birman’s Eqs. (8) and (9) for the case 
when JT = Ko. 

Then, since 

OT/Op = T/p 
Eq. (8) reduces to 
Ap/AW = p/W 
and Eq. (9) reduces to 
AV/AW =0 


These equations now are equivalent to those in the subject 
paper and are therefore not the result of neglecting the thrust- 
drag equilibrium requirement as claimed by Mrs. Birman but 
rather the result of an implicit assumption for the thrust varia- 
tion of a turbojet above 35,000 ft. Furthermore, despite Mrs. 
Birman’s having demonstrated one example where it is not strictly 
valid, this approximation still has justification in my mind for 
the purposes of the analysis in the paper. 
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Forced Convection Through a Laminar 
Boundary Layer over an Arbitrary Surface 
with an Arbitrary Temperature Variation* 


Myron Tribust and John Kleint 

University of California, Los Angeles, and University of Michigar 
Ann Arbor, Respectively 

August 11, 1954 


N A PREVIOUS report the authors considered the general prob- 

lem of forced convection in several systems when the wall tem- 
perature is nonuniform.! Table 1 summarizes the results of that 
report and is reproduced here because several new entries have 
been made in the table and a number of typographical errors 
have been discovered in reference 1. In all cases, if the wall 
temperature is given, the heat flux is to be calculated from 


x 
g(x) = h(é, x)dTw(£) 1) 
J é=0 
where 
q(x) = heat flux from the wall, B.t.u./hr. ft.? 
h(é, x) an integrating kernel (see Table 1), B.t.u./hr. ft 
F. 

7\,(x) = walltemperature, °F. 
x = distance along the wall, ft. 
E = dummy variable, ft 


and if the heat flux is given, the wall temperature is calculated 


from 
x 
TAs — ] = | o(&, x)g(&) dé 2 
J t=0 
where 
1 = reference temperature of fluid, °F. 
g(t, x) = integrating kernel (see Table 1), °F./(B.t.u./h1 


ft.*) ft 
Integral 1 is interpreted in the Stieltjes sense. 

In this note we give explicit directions for using the fourth 
entry of Table 1 to compute convection from an arbitrary sur 
face. The fluid properties are constant.** Aerodynamic heat- 
ing is zero. 

The method consists of first solving the momentum equation 
by using the velocity distributions of Hartree? and the ‘‘patch 
ing’’ technique of Eckert.* The energy equation is then solved 
using the method of Lighthill. 
the momentum equation are modified according to the method of 
Tifford® before being used in the Lighthill integrating kernel (line 
4, Table 1). The results of two typical calculations are shown in 
Figs. 1 and 2, where the experimental results of Giedt® and Drake 
method 


The shear stresses obtained from 


are given. We refer to the method as the H.E.L.T. 
(Hartree, Eckert, Lighthill, Tifford). A study of five references! 

will provide the background for following these instructions for 
making the calculations such as are demonstrated in Figs. 1 and 2 
For a given surface it is presumed the pressure distribution is 
known and either the heat flux or wall temperature are prescribed 
along the surface, with the unprescribed quantity to be deter- 


mined. 

(1) From the pressure distribution a calculation is made and a 
graph prepared showing u(x) and du/dx, where u(x) = velocity 
just outside the boundary layer in feet per second, and x = dis- 


tance along the surface in feet. 


* This work was done as a part of the Icing Research Program of the Engi 
neering Research Institute at the University of Michigan and was sponsored 
by the Aeronautical Research Laboratory of the Wright Aeronautical De 
velopment Center. 

+ Associate Professor of Engineering 

t Research Assistant. 

** See reference 18 for extension to air with variable properties 
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TABLE I 
INTEGRATING KERNELS FOR NONISOTHERMAL CONVECTION 
x 
ats) = (ME, a T9(€) TWw(t)—To"/ GE, x) (EE 
+0 
AUTHOR Er... SYSTEM __ME1 rHOD OF SOLUTION h(é,x) ake _gl€,x) — 
| Laminar flow over a flat plate | | velocity and temperature profiles | | 
| Zero pressure gradient. Fluid postulated linear in y. Thermal pA’? rey? ( it 3va)-1/3 | V3 -V2 | Mal-2/3 
RUBESIN 8 properties constant boundary layer thickness prepertional| i 7 | |eW7SMa7S Oscar Pr Rey |i- ny } 
| to momentum thickness. (Integral) | 
+ ——s inate naar amelie = + a 
| | Laminar flow over a flot plate Velocity and tenpercture profiles postu | | 
| Zero pressure gradient. Fluid | lated asa cubiciny integral method -" M4l-3 i v3 -2 3/4)-2/3 
ECKERT ® | Soeeuee anes aa 2338 pre [ifs] |awarteran o's Re,” | 64 
} | 
+ 7 . aie tt + _ — 
j | Laminar flow over a flat plate | Velocity profile taken as u*cy, not = | 
| Shear at the wall postuloted dependent upon x. Term vaT/dy v3i T 2.00 v3, -v3f 4/3 273 
LEVEQUE 10 | constont Fluid properties | Saaety Gogo ant of enarey eqseton | Stizgy: (P79) |(du/dy) 20) (x-€) Kear or p/%) [160/49] (re) 
conston resulting equation is solved | 
HH 4 — ——— — $$ 
Lominar flow over a surface with | Velocity profile token os uly | Bs en =23 
known variation in surface shear | where t(x)* wall sheor | Pi ek) ¢ { i - J Ig 
LIGHTHILL 4 stress. Constant fluid properties! Resulting differentia! equations | ad (rapt sid Lf ray] a2 \pm) [st Wr(z) ae | 
solved | 
T — rv ac: te 
“Wedge-flows’ Fluid properties | Velocity near wall taken os linear in y | tem¥"b— 5” D *Rey? [i-&F] | ms 
BOND Tr constant. (Velocity over wedge) Resulting differential equotion solved er ‘ ¢ yi ra?(s 2 [age] | 
given by u*qx™) | for two cases: (a) Step function in P be Sito) ce tiem leven em, x ) a; 
rempercture function “01 is 0 dmensioniess tabulated | - 
i pe Se 04 Step —— as a function of min reference 2. + a — 
| Laminar flow over a surface. | Velocity profile token os v3 | 
MODIFIED Fluid properties constant dr/dx small enough | els, J” 
LE VEQUE | | | de/dx smaii to fake vaT/ay negligible || Gs ne Te | 
incsdiioinsedl sheor 
+ ——T — — T + 4 
| Turbulent flow over a flat | Velocity and temperature - | 
RUBESIN } 12 | plote. Fluid properties con- profiles taken as following 2.0288k py v3 reine} » | _(2evi9s prs go? 08, zB (8; -& 
| stont V7 power tow. integral method | [32739)'(77ssoozesk ex * | 
— - _ ——— — — ena cee ee ee ee eae emanates _ — ee —| 
| Turbulent flow over a flat | Velocity profile token as /7 power | - 
SEBAN 13 plate. Constant fluid prop- | of y. Temperature profile token os 2.0289h py v9 reo? |. ary” | (8/90) Pr Ree? ,0% & hes | 
erties linear in y near the wall, \/7 power sa . | (e79nti79) oozeex | 
| | 
} >} —_—______— ——=} = — ene — 
| Experimental measurements on | Empirical equation “best-tit” to 
MAISEL AND 14 | @ mass transfer apparatus. | dato. (Transiated to air heat || oossk 08 [,_ ,gospou } 
SHERWOOD | transfer system by present x L a | 
| author) | 
+ + —__—— <1. - ie  iapaaperes — ————————+ - i 
| Laminar flow in o tube. Poro- | Differential equation solved by || * a ‘ eoBla-€ €) a ee. gs FY cemAl-€) 
GRAE TZ | bolic velocity distribution | separation of variables. Only | 3 a, TH3 4449 1138 a 
15 | Constant fluid properties first three eigen values known sii c one 0538 0488 | pte om ba ati 
2wep | Gn* 000 2564 8462 
4 + + -_ + : } . — —}— oie ~Alinn—e ——— 
| Turbulent flow of © liquid metal; Thermal conductivity of metal ifs vaP>') & 
POPPENDIEK | |6 | ina tube. Velocity profile postulated large enough to Sl2(P+i > e Pai F “| (x-€) 
| established Constant fluid render eddy heot transport neg- Psa 523)! a “g n Constonts cr, ond 
| t : we ond | 
| ee i. BE hn | age: Verey pome % wcsitas ws. UE IP tabulated in reference iz sua} ae ee 
POPPENDIEK Turbulent flow of o liquid metal | Thermal conductivity of metal poe! Pez - pi i : “95 ~pe 
| in @ tube. Velocity profile tulated large enough to render || __k — grr? (P+3)(Poi) = -~gyrr2 
MODIFICATION} 17 | established Gonstont fluid eddy heat tronsport negligible. (E*3)1 (P+2)* so (x-€) (Pe2) a EA I] Pe2ksPo (x-€) 
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Comparison of calculated and measured heat transfer FIG 


from a uniformly heated circular cylinder in cross flow.® 


& DISTANCE ALONG SURFACE OF CYLINDER 


2. Comparison of calculated and measured’ heat transfer 


from a uniformly heated elliptical cylinder in cross flow. 
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Fic. 3. Function to be used in calculating with Eq. (2 


2) Starting from the stagnation point, it will be found that 
for a certain distance the graph of u(x) versus x is a straight linc 
In this region the momentum thickness J, in feet, is given by 
0.2921x (Re,) 


formed from the local velocity, U(x 


i = where the Reynolds modulus, ke 
, and the distance x from th 
“wedge parameter,” #7 


stagnation point In this region the 


1 (see reference 2 ) 
(3) Downstream of this region it is necessary to solve the fol 


lowing equations (see reference 3 


di? l—m # du 
dx m u dx 
#(du/dx 2m 
= F(m)|? (2 
v 1+ m 


F(m 


Thus, when the velocity is 


kinematic viscosity in feet per second function 
obtained from references 2 and 38 
no longer linear in x, we use the known values of J?, du/dx, to 
calculate F(m) and from the graph, Fig. 3, we find m 
Eq. (1 
downstream 
F(m), ete 


From 


then, a new value of J? is computed for a position Ax 


In this way a tabular set of values of 3(x), m(x), 
is prepared. Eckert* suggests the isocline method 
of solution 

(4) The wall from the equation 


shear stress is calculated 


T u(du dy)y 0 MUX F(m)f"(O, m)/d(x 


Fig. 4 gives a graph of f"(O, 17 

(5) Before substituting this shear stress in the Lighthill inte 
grating kernels (line 4, Table 1), the correction of Tifford is 
applied, which may be written: 


} 2m (Fm 
Pr 


Py = Prandtl modulus of the fluid, dimensionless 
This is a quasi-empirical correction to Lighthill’s approxima 
tion for the more complicated velocity distribution in the bound 


ary laver 


TARLE 2 


The Functions F(m) and f"(0, m 


m F(a1)}? f"(O, m) 
0.0 0.4696 0.470 
0.1 0.410 0.674 
0.2 805 
0.3 0.900 
0.4 0.977 
0.5 1.040 
0.6 1.095 
0.7 1.134 
0.8 Lz 

0.9 1.20 

1.0 1.23 
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(6) Values of 7» (effective) are then substituted into the int 
grating kernels of Table 1, line 4, and the integrations performed 











graphically Table 2 gives numerical values for F(m), f’(0, 
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inte 
cined } Ducts with Noncircular Cross Sections 
E.R. G. Eckert* and T. F. Irvine, Jr.t 
Transfer Laboratory, University of Minnesota, Minneapo 
S August 12, 1954 
Ann 
1 } Des HeAT TRANSFER Laboratory at the University of Miu 
‘ nesota is engaged in an experimental program to study the 
flow characteristics in noncircular ducts Flow visualization 
La techniques and probing of the velocity field have been used to 
= examine the flow of air through a triangular duct having the shap: 
of an isosceles triangle with a height to base ratio of 5:1 Phe 
>) duct length was 80 hydraulic diameters and the flow entered the 
duct through a well-rounded entrance 
' hese studies revealed the unexpected fact that over th 
283 whole tube length laminar and turbulent flow coexist side by side 
over a considerable Reynolds Number range The flow near the 
] . 
we ' b f the tube is found to be turbulent, whereas near the apex Fic. 1. Visualization of flow through triangular duct by 
Atica F siabiiiialnd Meinl smoke. Flow from right to left Smoke is ejected through small 
ow remains la lay stb : “sa : 
: a probe lop picture indicates laminar flow Center picture 
ts on lo investigate these phenomena the flow was made visible by begianing instability Bottom picture turbulent flow 
>». 41, injection of smoke through a small probe which could be moved 
sie } along the height of the triangular cross section As the probe 
was moved from the apex toward the base it was first observed that 
nsfe ; T tion to turbulence in triangular tube 
“bie over a certain distance the smoke appeared as a completely stag LO ransition TO Ti = 
y of nant thread. As the motion toward the base continued the smok« . “7 
apr 7 si Flow Cross-Section 
thread began to exhibit sinusoidal waves (similar to Tollmien <— 
raw Schlichting waves traveli ig dow nstre um und fluctuating within ~ogl i - Re=- Vi Dh 4 
the lengthwise plane which contains the height of the cross sec < vy 
i, R : . 3 ® x 
tion. Finally, as the probe reached a position closer to the base, o. V; =Mean Flow 
the amplitude of the waves increased until at a certain int the So L L 4 
Uni ie : a ¢ eed r — : i che i po rs Velocity 
smoke thread dillused within a short distance downstream from : 
ty of aap S D,=Hydraulic 
j the probe, indicating fully developed turbulence ~ 
: =o4l Turbulent Diameter j 
ure Fig. 1 consisting of three separate photographs shows th« Z d 4 
on ‘ +e ® = 
smoke thread as observed at three locations within the tubc In oO one vy =Kinematic 
, this figure the base of the triangle is at the bottom of the phot ) S02 Viscosity 
orf 2 L , 
in graph. In the top picture, the completely straight smoke trail <= Laminar 
} is indicative of laminar flow, the center smoke trace, which is neat oO Zone idddddaccacaacccdarr 
bu the edge of the turbulent region. is showing the long wave-length 1 1 A , ne 
disturbances mentioned above, and the lower trace is rapidly O 2000 4000 6000 8000 
dispersed by the turbulent flow existing near the triangle base Re 
ance From such flow observations the extent of the laminar regim« Fic. 2. Laminar and turbulent flow regimes in triangular duct 
ries ° 
could be determined is determined by flow visualization 
Fig. 2 presents the results of a large number of such observa 
ee tions in the 5:1 duct It is interesting to note that the Reynolds 
) Number for which the flow is laminar in the whole cross sectio. is f f 
dary quite low in comparison to values for a round tube and that, on 14 Velocity profile a ong center line 0 duct , 
pod : the other hand, the laminar range which decreases in size with | Flow Cross-Section 
increasing Reynolds Number can still be observed for quite high | | 
i Seta ‘big ae | er - 12°— / | 
Reynolds Numbers. In this figure the hatched section indicates | x 
the spread of the data obtained in the above described way — | 
“| | naan ape ig 1.OF L : 
| rhe velocity field in a cross section 80 hydraulic diameters 
} downstream from the entrance has been measured by pitot 08! Re=6,600 “ 
- . 5 ° = | 
i tube traverses for various Reynolds Numbers V a ° Y= Point Velocity 
i Fig 3 shows as an example the velocity distribution along the V; O6| ig V;= Mean Flow ! 
height of the duct from base to apex l Velocity 
i] 
In the laminar range near the apex the velocity increases with 0.4} Vi ViDs 
the square of the distance from the apex in agreement with cal | 
culations for laminar flow near a corner. The extent of the 
Pa 
: : ‘ = | 
laminar regime as obtained from smoke observation ends at a 0.2 Q Dy, Hydraulic D ameter 
location indicated in the figure as point ‘‘a.’’ The part of the v= Kinematic Viscosity 
: ” Coes ° 4 rn i 
profile between ‘‘a’’ and ‘‘b’’ continues, however, to follow a O 02 04 0.6 08 LO 
parabolic variation. A critical Reynolds Number based on the . X 
rng : ak eee Distance from apex “/, 
velocity averaged over the duct width at location ‘‘x,’’ and on 
“To channel width at this location has at point ‘‘a’’ the value 102 and Fic. 3. Measured velocities in triangular duct. Point “a’ 
location of beginning instability (center picture of Fig. 1 
_ i * Professor of Mechanical Engineering Point ‘‘b’’—flow starts to deviate from calculated laminar pro 
T Instructor, Mechanical Engineering Department file 
} 
b 
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at point ‘‘b’’ 600. It remains to be investigated how far these 
values are dependent on the apex angle of the duct cross section. 

The program” is being extended to study the heat-transfer char- 
A detailed investigation of the turbu- 
The purpose of this note is 


acteristics of such ducts. 
lence characteristics is also planned. 
to report on the coexistence of the two flow regimes described 
above. 

* Sponsoring agency, Wright Air Development Center, Contract No. AF 
33 (616) 474. 


Hypersonic Similarity and the Tangent-Cone 
Approximation for Unyawed Bodies of 
Revolution; 


Ronald F. Probstein and Kenneth N. C. Bray 

Assistant Professor of Engineering and Applied Mathematics, 
Brown University, Providence, R.I., and Assistant in Research, 
Department of Aeronautical Engineering, Princeton University, 
Princeton, N.J., Respectively 


August 18, 1954 


_ IDEA OF THE tangent-cone approximation, as applied to 
determining the pressure distribution on a pointed unyawed 
body of revolution in supersonic flow, is to replace each trans 
verse section of the body by a section of a coaxial cone tangent 
to the body surface. The pressure is then assumed to be that 
which would exist on the cone, at the same free-stream Mach 
Number. Althougi: difficult to justify rigorously, this theory 
has been shown by Halbmillion and Kulishek,! among others, to 
give good results for slender ogives when compared with exact 
rotational characteristics calculations. All previous workers 
utilizing the tangent-cone approximation have determined the 
cone surface pressure from the numerical solution of the Taylor- 
Maccoll equations for supersonic conical flow (e.g., reference 2), 
and were therefore unable to obtain an analytic representation of 
the pressure distribution 

According to hypersonic similarity,* * the pressure on the sur- 
face of a series of slender unyawed bodies of revolution, with the 
same thickness distribution along their axes, may be written 


p =e 

— = G(K, £) (1) 

pi 
Here, K = M,/(l/d) is the hypersonic similarity parameter, / is 
body length, d is body diameter, J; is the undisturbed free 
stream Mach Number, and & = x/l, where x is axial distance 


measured from the nose. Also, the drag coefficient based on 


maximum frontal area takes the form 
D F(K) 


Cp = a ma yA 
5 pM2(2 a?) Mi? 


Hypersonic similarity has been found to apply not only for 
high Mach Numbers, but also at moderate supersonic speeds,*: 
Lees? making use of the hypersonic similarity 


6 


(M, = 2.5 — 3). 
concept has determined an approximate expression for calculating 
the pressure on the surface of a slender unyawed cone, when the 
conical shock is not too far from the surface. The present note 
will show by the use of Lees’ equation and the tangent-cone ap- 
proximation, that a simple and accurate relation, which depends 
only on the local body shape, is given for the pressure distribu 
tion on a sharp nose unyawed body of revolution in supersonic 

t The present work was sponsored by the Aeronautical Research Labora 
tory, Wright Air Development Center, United States Air Force, under 
Contract Number AF 33(038)-250, and was completed while the senior 


author was at Princeton University. 
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flow when K > 1. 
pressure distribution can be integrated analytically (or numeri 
cally if need be), and the overall wave drag determined. 

It may be argued that for values of the hypersonic similarity 


parameter in excess of one, Eggers et al.7~* have shown that the 


surface pressure and Mach Number distribution about an un 
vawed body of revolution can be calculated by assuming a conical 
shock at the nose followed by a Prandtl-Meyer expansion. Al 
though this method gives good results, the algebra is complicated, 
as conditions at a point on the body do not depend solely on the 
More recently, Van Dyke’ has 
theory for slender 


local body shape at that point. 
developed a “‘hypersonic small-disturbance’’ 
three-dimensional bodies. The present authors feel that the 
complexity of this analysis make its use in preference to the 
method of characteristics questionable, particularly in view of the 
restricted accuracy 

For slender bodies and the lower supersonic speeds the linear 
ized Karman-Moore!! theory offers a good approximation to the 
pressure distribution on unyawed bodies of revolution. Van 
Dyke’s!? second-order theory increases its accuracy in the range 
where the hypersonic similarity parameter is less than unity 
However, the algebra is somewhat cumbersome, and in general 
the results cannot be expressed in an analytic form. In the 
present note it will be shown for this case, when the hypersonic 
similarity parameter is less than unity, and the shock or Mach 
cone is far from the surface, that application of the tangent-cone 
approximation utilizing the Karman'* slender-body linearized 
cone solution yields a simple and accurate analytic representation 
of the pressure distribution and drag 

Lees’ showed for slender cones when the conical shock is not too 
far from the cone surface, and when hypersonic similarity applies, 


that the surface pressure is given by 


on 
__ 7 ee —— (K,* — 1) + 
pi yrs 
+ ] 
7¥ (Ks — K 2d : i J 
ly —14+(2 K,2)4 
where K, = 16, (6, is the cone half angle), and 


7¥ + !1 a 2 9 
K, = M0, = - = K. + (: ) x2 4. = { 
.. = \ ¥+3 Y¥+3 


(@; is the conical shock half angle). These relations were obtained 
by considering the limiting forms of the Taylor series expansions 
for the solution of the differential equation for conical flow behind 


the shock, the boundary condition at the shock surface, and the 


oblique shock relations. The approximate limits of applic 
ability of these relations is given as @, < 20° for air, 3 < /, - 
and K- > 1. 

In the tangent-cone analysis, 0, of Eqs. (3) and (4) is replaced by 
dy/dx, the local slope of the body, so that Ke = W/,(dy/dx 
For a family of similar bodies of revolution with different fineness 
ratios, K, = A,(K, &), in which case from Eqs. (3) and (4) it is 


clear that similarity as expressed by Eq. (1) is satisfied and the 
problem is determined. 

By way of example, pressure distributions over circular-ar¢ 
ogives will be determined, because the rotational characteristics 
solutions are available for comparison with the approximate 
formulations to be presented. Now ogives of different fineness 


ratios are not strictly similar bodies. However, hypersonic 
similarity requires fineness ratios sufficiently high, so that within 
the limits of applicability of the hypersonic similarity law the 
ogive approximates a similar body very closely (this is borne out 
by the results of reference 6). Thus, for a slender ogive with 


(1/d)? > 1 it is not difficult to show that 


+ 
} Apart from what appears to be an incorrect scale in the comparison curv 


of Fig. 7 in reference 10, it seems to the present authors that there would be 


an even larger error than indicated This stems from the fact that at the 
leading edge the pressure coefficient cannot be the exact value as is shown 


by the results of Fig. 4 in the same reference 


Furthermore, once the shape is specified the 
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K K (1 — &) (5) 
[his is the equation for a parabolic body of revolution, and for 
values of the fineness ratio greater than 3 (approximately ) may be 
taken to represent the equation of an ogive. It is to be noted 
that the minimum fineness ratio of 3 corresponds to a maximum 
nose angle of about 20°, which is the limit given by Lees in con- 
nection with his hypersonic similarity analysis. Of course, for an 
ogive, Ke ~ 0 as § > 1, so that Lees’ assumption of K, > 1 can 


never be satisfied near the rear of such a body 
The wave drag coefficient for a body of revolution is 


Ps 
Cp = 2 | Cpn dn (6 


0 
where » = y/(d/2), and Cp is the pressure coefficient. However 
for an ogive, from Eq. (5), n = #2 — &), so that 
Py 
Co = 4 Cypt(1 — §)(2 — &) dé (7) 
0 


To simplify the writing y is taken equal to 1.4, so that the pres- 


sure coefficient from Eqs. (3) to (5) is given by 


=(? 2 ) se 


y 


CpM\? = 1.0415K%(1 


— £)? — 0.4545 4 


1.6563AK%1— &)? (8) 


V 1.0842K41 — &)4 4 
If the above relation is substituted into Eq. (7) and the integra- 
tion carried out, the following functional relationship for F(K), 
which satisfies the condition of hypersonic similarity [Eq. (2)], is 
obtained: 
= CpM,? = 0.3471K? — 0.4545 — 1/K*}4.2780x3 — 


0.4641) In (2.3998x 4 


F(K 


1.3816K2 + 1.0554)?x + (0.6076K? 4 


1.3089K2 + 1)} (9) 
where x V 0.2975K4 + 0.4545K2. Again, this formulation is 
expected to be valid for K > 1, with (1/d 3 and M, > 3 (ap 


proximately 
At very high Mach Numbers, the value of the pressure coef- 


ficient on a slender unyawed cone is given by (see, e.g., reference 


WU 5*Co.... = : = K.? = 2.08K,* for y = 1.4 


rhis is seen to be close to the y = | Newtonian value of 2 K,? 


From Eqs. (5) and (7) the drag function for air becomes 
208 
CpM 2 = “ K? 10 


o 


F(K) = 


The above result corresponds to the Newtonian value in which 
the influence of finite surface curvature inducing a centrifugal 


effect is ignored Note, of course, that the foregoing relation 


uuld have been obtained directly from Eq. (9) as the limiting 
case where A? >1. One might therefore expect for a value of K 
greater than about 3, that Eq. (9) could be replaced by the New 
tonian value with only a small error 

When the hypersonic similarity parameter is less than unity, it 
is clear that the previous formulation is invalid 
body is slender and the Mach Number moderate, then the shock 
or Mach cone is located far from the body surface, and it should 


the 


How ever, if the 


be possible to apply the tangent-cone approximation to 


Karman slender-body result for the linearized supersonic flow over 


cone. As given by Karman!’ 
Pe a 2 a 2 
— 1= ¥K,? In - = yK,? In — 
Pi 6. M,? — 1 Ke 
rherefore, from Eq. (5) the pressure coefficient on an ogive is 
» ») 
2/(?P oa od - 
( — | = C, Vi: 2A%1 — &)? in 1] 
y\p~r A(1l —¢& 
his leads to the following equation for the drag function 
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Fic. 1 Pressure distribution on ogives by tangent-cone Lees 


approximation and rotational characteristics for values of the 


hypersonic similarity parameter K = 1.0, 1.5, 2.0, and 2.29 
2 ys a) 
F(K) = CpM? = ~K? In— 4 K? 12) 
3 K 18 


In Fig. 1 the pressure distributions for the range where A > 1, 


calculated by means of Eq. (8), are compared with those ob 


tained from rotational characteristics solutions, for values of AK = 
1.0, 1.5, 2.0, and 2.29. All 

taken from Rossow® with the exception of the A = 2.29 calcula 
tion (ref. 14, M@, = 8, 1/d = 3.5). 


hypersonic similarity to the characteristics 


the characteristics solutions were 
The @ priori applicability of 


solutions has been 


assumed, in that only the faired values from Rossow, which en- 
compass the Mach Number range of 3 to 12 and fineness ratio 
range of 3 to 12 are shown It can be seen, as might have been 
expected, that the agreement is good except over a small portion 
2 presents the rotational charac 
K = 0.5, 0.75, 


comparison with those calculated by Eq. (11 


of the body near the rear. Fig 
and 1.0 in 


Away from the 


teristics pressure distributions for 


1, the agreement is again seen to be good 
CpM,? 


> 


is shown plotted in Fig. 3 


rear, and for K - 
The drag function F(K) = (10), 


are six points 


is given by Eqs. (9), 
and (12 Also shown 
representing the drag parameter computed by the use of rota 
tional characteristics theory All the values, with the exception 
of the K = 2.29 point, which was computed using the pressure 
distribution of Fig. 1 and Eq. (7), were taken directly from Fig 
this case the minimum VW, is 3, 
while the minimum value of //d is also3. What is clear from this 


figure is that for K > 1 the expression of Lees in combination with 


7 of reference 6. In value of 


the tangent-cone analysis gives an excellent approximation to the 
drag coefficient, while below a K of one the Karmén tangent 
cone relation is excellent. If both formulations are taken to- 
gether one has an analytic representation of the pressure dis- 
tribution and the drag coefficient for unyawed slender pointed 
bodies of revolution over practically the entire range of super- 
Of course, this has only been checked for 
Mach 


sonic Mach Numbers 


ogives with fineness ratios greater than 3 and free-stream 
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Fic. 2. Pressure distribution on ogives by tangent-cone 


Karman approximation and rotational characteristics for values 
of the hypersonic similarity parameter A = 1.0, 0.75,and 0.5 


Numbers greater than 3, but it can be expected to hold for other 


body shapes as well 
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An Application of the Buffon Needle Problem 
to the Determination of the Probability of 
Obtaining a Free Flight Arrested Landing 
in Carrier Landings 


James E. Martin and Walter W. Hoy 
Chance Vought Aircraft, Incorporated, Dallas, Tex 
August 26 1954 


REE FLIGHT ARRESTED LANDINGS of carrier based aircraft are 

those in which the arresting hook engages the cable before 
either the main or nose gears contact the deck In designing an 
arresting hook for a particular aircraft configuration, it is desir 
able to know the probability of the occurrence of such a landing 
condition. This problem has been solved by application of the 
Buffon Needle problem. 

The Buffon Needle problem can be stated as follows: What is 
the probability of a needle, of length f, intersecting one of an in 
finite number of equidistant parallel lines when dropped from an 
irbitrary height above the lines? It is noted that the length 
is considered to be less than the width of strip between two ad 
jacent lines 

To adapt the Buffon Needle problem to the problem at hand, 
consideration must first be given to the “length of the needle, 
or the effective hook length. This length will vary depending 
on the airplane’s sink speed, engaging velocity and attitude, 
and the ‘‘needle rotation”’ will be determined by the airplane’s 
vaw angle. The set of parallel lines will of course be the arresting 
cables aboard the carrier with due consideration to there being 
only a finite number of cables 

In order to evaluate the effective hook length and probabili 
ties with respect to a carrier landing, certain variables, and their 


1 


distributions must be considered. The notation and parameters 


to be considered are as follows: 


t engaging velocity, ft. per sec 

Z sink speed, ft. per sec. 

A attitude, deg 

20 height of main gear above the deck when the hook contacts 
the deck and the airplane has attitude @ 

y yaw angle, deg 

h effective hook length, ft. (i.e., the distance the hook head is 
on deck, measured from its deck contact point to a point 
corresponding to the time the miain, or nose, gear makes 
deck contact 


h > [x z cot 30 4 
f, (0) frequency function of @ 
fo(z) frequency function of z 
(x) frequency function of x 
FF free flight arrested landing 
P(FF\6) = the probability of a free flight arrested landing given an 
attitude 0 60 
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Phe variables «, 2, and @ are assumed independent, in a probability 


ense i.c., the joint frequency ¢g(i, 2,6) = fifof The analvtical 
expressions of the functions f; may not be known; however a good 
ipproximation may be made using test data and a sample rela 
ca may be found analytically 


tive frequency plot The variable 


is a function of #, 2 and @ or may be empirically determined by 
measurement from scaled drawings or prints 

Making use of the Buffon Needle problem, the probability of a 
free flight arrested landing, for a particular flight condition and 


n infinite number of cables, may be evaluated and has the form 


J @ sin y¥ 
P, = — cot (30° — @ 
id d y 


Since only a finite number of cables exists, the probability, p, of 
landing in the area of deck containing cables must be considered 
Thus the probability of a free flight arrested landing aboard a 
carrier becomes ?, = pP,. 

The probability of this type of landing for a particular air 
plane may now be determined. With the assumption of inde 


pendence and knowing /i(@), f2(2), and f;(@), this probability takes 
the form 


. . 
P, = { | { fi (0) f(2)f3(4)P(0, 2, «) dd dz di, 
JRiIR SIR 


} 


where R,, Rs, and R; are the ranges of definition for each of the 


variables. In general, ?; must be evaluated by considering the 
definition of an integral, that is, a limit of a triple sum. It is 
noted that the accuracy desired is dependent upon the limits of 
summation. The use of IBM equipment will facilitate this cal- 
culation 

The results of a study of this type may best be presented 
graphical The probability P2 may be plotted as a family of 
curves depending on #/2 and @, thus providing a rapid means of 
obtaining 7. for a particular flight condition. By summing P 
first with respect to 2 and #, the P( FF\@)) may be plotted as a 
function of @, and similarly for P(FF\z) and P(FF\z). Final 
variable provides 


summation on the remaining independent 


the value of ?; for a particular airplane. 


Parameters Predicting the Initial and Final 
Collapse Pressures of Uniformly Loaded 
Spherical Shells 


Bertram Klein 
Research Engineer, University of California at Los Angeles 
September 2, 1954 


HE PROBLEM OF THE COLLAPSE of a spherical dome under 
uniform external pressure was considered solved! until re- 
cently when tests were undertaken to determine the collapse 
pressures of fairly flat shells.2 There was a rather large scatter 
when the results were plotted against the apparently correct 
parameters defining the problem. It is the purpose of this note 
to show how the choice of more rational parameters reduces the 
seatter considerably and at the same time allows the designer 
to account for initial manufacturing eccentricities in choosing 
the allowable collapse pressures 
As already established,' one parameter is R/t, which defines 
the flexibility of the shell in resisting external pressure. The 
other parameter is related to initial weak spots in the structure 
which always exist in practice. It is well known that the pres 
ence of these initial irregularities is an important factor in the 
collapse strength of thin shells.*~* In the present instance, it 
is believed that both thickness and curvature play an equal role 
in determining the importance of initial irregularities. There- 


fore it is suggested that e/V/ Rt be used as the eccentricity param- 
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Fic. | Initial collapse stress of spherical shells under uniform 
external pressure 
eter. Here ¢ is taken as the average maximum initial departure 


from the mean radius of the shell 
In the past it has been considered natural to plot A vs. R/1, 


where A is the variable coefficient in the equation 
o0 = KEt/R 


However, it has been found that using A(t/R o/E, reduces 


the test scatter. Figs. 1 and 2 show the variation of the initial 
and the final collapse stresses with the chosen parameters (plotted 


vs. ¢/E,) for known test data Where R was not specifically 


given, it was computed from the relation: R = (8/h) + (h/2), 
where & is the base length of the spherical segment, and h its 
height Notice that the angle subtended by the segment is not 
considered to be a parameter in the problem. This is true as 
long as its value does not become too small 

The shapes of the curves in Fig. | are similar to the one shown 
for a perfect shell (A = 0.606 Also it is seen that the drop in 


initial collapse stress with e/ VY Rt is very rapid at first, then 


tapers off. From Fig. 2, it appears eccentricity is not significant 
in determining final collapse pressure unless the eccentricity is 
sizable 

No doubt considerable additional data are available which can 
be used to establish more securely the authenticity of the curves 
The author would be interested in learning of such data and how 


it fits the plots 
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Boundary-Layer Stabilization by Surface 
Cooling in Supersonic Flow 


E.R. Van Driest and J. Christopher Boison 

Aerophysics Department, North American Aviation, Inc., Downey, 
Calif. 

August 24, 1954 


A PROGRAM* IS BEING CONDUCTED in the Aerophysics Depart 

ment of North American Aviation, Inc., to study experi 
mentally the effect of surface cooling upon the stabilization of 
boundary under wall-to-local-free 


supersonic laminar layers 


stream temperature ratios normally encountered only in transient 
flight. 
A 20 


duction of gaseous nitrogen, the prime coolant being liquid nitro- 


(included angle) 6-in. cone is cooled internally by intro 
gen. By this method, steady surface temperatures at any level 


down to about —300°F. can be attained. The studies are being 
carried out in the NAA 3! 
sonic wind tunnel operating in a Nh 


3.70 at Reynolds Numbers up to 10° per inch with 70°F 


by 3'/.-in. continuous-flow super 
fach Number range of 2.45 to 
supply 
temperature 
Figs. 1 and 2 are typical schlieren photographs for a smooth 
cone. Fig. 1 shows the cone, uncooled, with transition of the 
boundary layer about half-way down the cone. The position 
of the transition region was also verified by surface-temperature 
measurement. The Reynolds Number of the entire cone was 
54 & 10°. The local Mach Number just outside the boundary 
laver was 3.30. Fig. 2 shows the cone under the same flow con- 
ditions, but cooled internally to approximately —200°F., corre 
sponding to a wall-to-local-free-stream temperature ratio of 1.5 
This temperature ratio is theoretically sufficient to stabilize the 
* This research was supported by the United States Air lorce, through the 
Office of Scientific Research of the Air Research and Development Com- 


mand 
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Re cone *9.4x108 


Fic. 1. 








ns COOLED CONE SOUNDARY LAYER. 
Fic. 2 
laminar boundary laver for all Reynolds Numbers A lami 


nar boundary layer over the entire length of cone is indeed evi 
dent in Fig. 2, thus indicating a delay in transition due to cool 
ing. It is interesting to note that in this particular set of tests 
in which a reflected wave intersects the aft end of the cone, the 
the intersection, when 


boundary layer remains laminar, after 


the surface is cooled Although the Reynolds Numbers of the 


cone are not large, the experiments can be quite instructive, 
especially in the study of roughness effects 
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Errata—‘‘Simplified Analysis of Turbulent 
Boundary-Layer Development Along 
Cylindrical Surfaces with Variable Free- 
Stream Mach Number’’' 


Hans U. Eckert 

Air Research and Development Command, Wright Air Development 
Center, Dayton, Ohio 

October 29, 1954 


The fifth word 


P** 696, column 1, line 11 from the bottom: 
The third 


should Column 2, Eq. (9 
term should have a rectangular bracket on the right side preced- 


read ‘“‘theorem.”’ 
ing the equals sign. 
Page 698, Eq. (22): The last factor should have a minus sign 


in the exponent. Column 2, line 13: ‘1/6"’ should be inserted 
Replace ¢ with p 


Column 2, line 16 


before the integral. Eq. (27 
Page 699, Eq. (28 Replace ¢ with p 
from the bottom: The term in the center should read cy/cy; 
Page 702, Eq. (43): Both integrals of the function f are part of 
the exponents and should be raised slightly 
Page 703, column 1, line 10: A should be a lower case letter; 


kandk + 


1 are subscripts to x. 


Page 794, Eq. (27a), column 2, top line: Replace ¢ with p 
Page 705, column 1, line 11: Instead of ‘Eq. (2.4a)"’ substi 
tute “‘Eq. (2.la).” Line 24: Instead of “Eq. (2.1 substi 
tute “Eq. (2 4).” 
REFERENCE 
' Eckert, Hans U Simplified Analysi Turbulent Boundary-Layer 
Development Along Cylindrical Surfaces with Variable Free-Stream Mach 


Vumber, Journal of the Aeronautical Sciences, Vol. 21, No. 10, pp. 695 
706, October, 1954 
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Beam Strength and Curvature Under Combined 
Tension and Bending in the Piastic Range 


Anthony J. Barrett 

Technical Assistant (Structures), Royal Aeronautical 
London, England 

September 9, 1954 


society, 


nana AND AND ROACH evaluate interaction curves between 

bending moment and axial load for a beam of rectangular 
cross section in the plastic range, making approximate allowance 
for section deformation The Ramberg-Osgood form of the 
stress-strain curve which is used is given in terms of true stress 
It is of interest to examine the errors involved if the effects of 
section deformation are ignored and the more convenient appar 
ent stress is used 

When bending moment is zero, the values of permissible axial 
load will obviously be as correct when evaluated as the apparent 
stress multiplied by the original cross-sectional area as they 
would be when evaluated as the true stress multiplied by the de 
This assumes, of course, that the allowance made 
When bend 


ing is present, together with axial load, no error is incurred by 


formed area 


for section deformation in the latter case is exact 


neglecting section deformation in the evaluation of the force con 
tribution of an elemental area within the beam; excepting insofar 
as the stress value used will be slightly incorrect, since the strain 


value at that element will have been evaluated on an erroneous 


distribution. The distribution constants are 
co = (1/2) (er + e2) + (1/16)ey(er — 
a. = (1/2) (¢ — @€o) + (1/S8)ey-(e1? — e2? 


ind if the effect of section deformation is neglected only the 
first term on the right-hand side is considered. In the extreme, 
when axial load is zero, ¢; —¢.. Ignoring section deformation 
therefore incurs a very small error in c; and an error of magni 
tude (1/4)eye,? in « The evaluation of the bending moment is 
further prejudiced by ignoring the movement in the centroid of 
the cross section The overall effect of these errors, in the ex 
treme case of pure bending, is best appreciated if the second terms 
in the expressions for cy and ¢; and the movement in the centroid 
are ignored and the bending moments recalculated for compari 
son with the values given by Frankland and Roach. Their 
integrals B, C, and D, simplify to the functions presented in refer 
ence 2 as $3, @) — di, and @ and the expressions for bending mo 
ment and axial load are also simplified to those presented in refer 
ence 2, after taking account of a variable maximum fiber stress 
and the different notations employed 

A numerical comparison is made in Table 1 to three significant 
figures. The modifications brought about by consideration of 
section deformation is only of the order of one per cent, even at 
As the basic 


the high plastic strains represented by m = 16 
assumptions of the method hardly justify the use of interaction 
curves to such precision, it is doubtful whether it is necessary to 
incorporate the effect of cross-section deformation for most prac 
tical purposes. The use of apparent stress and ignoring cross 
section deformation simplifies the calculations and should 
provide data as reliable as can be expected in view of the basic 
assumptions of the method 
A nondimensional curvature «x is evaluated! as: 


and this is a sufficiently accurate expression when the strains are 


small by comparison with unity. More nearly 


€ = 
kK * 
oT éy(¢ r ¢ 
lhe use of the approximate expression can involve errors in curva- 
ture of the order of 5 per cent when m = 16. Although in the 


practical ranges of maximum strain, the use of the approximate 


FORUM il 


raBLe | 
W/ My (for P/P 0 
Material Measure of Ignoring 
Characteristic, Plasticity From Curves of Cross-Section 
n ” Reference | Deformation 
1 1.0 1.40 1.40 
$4) 1.70 1.71 
16.0 1.9% 1.9S 
15 1.0 1.39 1.39 
1.0 1 64 1.64 
16.0 1.82 1.84 
21 1.0 1.38 1.38 
0) 1.59 1.60 
16 we 1 73 


expression for curvature is probably quite satisfactory, the errors 
involved are greater than those involved in neglecting the effects 
of section distortion. 

The restrictions on the use of the method for curved members 
under axial compressive load may be illustrated by the presen 
tation of curvature, strain, and end-load relationships in the form 
of Fig. 1, where d?w/dx? is the change in curvature from the un 
loaded form and 2D is the beam depth in the plane of bending 
The curves are appropriate to the standardized material pre 
viously used,” apply to apparent stress, and were calculated with 
out consideration of section deformation For the reasons dis 
cussed in reference 2 the maximum stress considered is the 0.5 
per cent offset yield stress (f// 1.05 It is convenient to 
specify an equivalent eccentricity of axial load, n, as the ratio 
of bending moment to axial load 

Assuming the centroidal axis under no load is given by 


ay cos (rx /1), and that under load it is given by a cos (74x 
then 

d*w /dx? r/l a= 
at the critical mid-length section Now a is the eccentricity 7 
and 

dw rD n a 

D 

dx? l D D 

which is the equation of a straight line such as XX on Fig. | Phe 


slope of the line is proportional to the elastic instability stress of 
an initially straight member of the same dimensions and the inter 
cept on the 7/D axis is the nondimensional initial bow a)/D. As 
the curvature increases the intersections of the line and the con 
tours of stress f4/fo and fp/fs record the stress levels in the inner 
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and outer fibers, respectively; f. being the 0.2 per cent offset 
vield stress 

For elastic strain reversals an equation of the Ramberg-Osgood 
type still gives an adequate representation of the relationship 
between stress and strain. The application of the method is 
therefore valid as long as any strain reversal which is taking place 
within the member is elastic. Referring to Fig. 1, the applica 





tion of the method to an initially curved member under axial 
load is valid for characteristic lines relating D(d?w/dx?) and n/D 
which do not cross stress contours, greater than that appropriate 
to the limit of elasticity, in successively decreasing magnitude Fig. 1 
For the material used in the calculation of Fig. 1, values of f/f» 
less than 0.75 are essentially elastic so that the use of the method 
for members in this material having characteristic lines passing where Y and JY are the forces in the directions of — and », respec- 
through the region shown hatched would not be valid Initially tively 
curved compression members, with initial bow of the order of Upon substitution for the pressure p in Eq. (2) from the mo 
beam depth 10 and greater, would have characteristic lines mentum equation in the form 
passing clear of the hatched region and this value gives a rough P tdW 
idea of the range of applicability of the method. “ os + P = const. 
ol 2d¢ dt p 
REFERENCES the equation for the forces on the cylinder becomes 
re N : ¢ Strength vile) ombined Tension , T’ 

and oe a cay dene gat dee of st pelletesoe Jbseiiat P ; me (d + 7dé ~ 1 + id¢) 
Vol. 21, No. 7, pp. 449-433, July, 1954 : é —— ao 

* Barrett Anthony J Unsymmetrical Bending and Bending Combined (4) 
vith Axial Loading of a Beam of Rectangular Cross Section into the Plastic 
Range, Journal of the Royal Aeronautical Society, Vol. LVII, pp. 503-509, = The velocity potential ¢ at the surface of the cylinder is from 
August, 1953 Eq. (1) 


‘ ro asin@é—y 
g = 2Uacos@ — tan ~! — 
2x | acos@é@—x 
asin@é + y asin 6 + (a*y/r?) 
tan! ; + tan™ 
acos@é—x a cos 6 — (a?*x/r?) 
The Unsteady Forces on a Circular Cylinder in ee sin 0 — (a2y/r2) ]) (f 
an r 


the Presence of Two Symmetrically Disposed a cos 0 — (a2x/r?) 
Line Vortices 


The velocity potential for a cylinder of varying radius is useful 
in slender-body analyses. Eq. (5) applies to such cases.  Be- 
S. Sherman Edwards : : 
Aeronautical Research Scientist, Ames Aeronautical Laboratory, ; 7 ee a ike ; : aad 

NACA. Moffett Field. Calif. tions, Eq. (5) with independent variables T, a, x, and y can be 
October 95 1954 differentiated with respect to time and substituted into Eq. (4). 

' 


cause of the image vortex representation of the boundary condi- 


Upon integration of the resulting expression around the boundary, 


name . . the force on the cylinder is 
ee has shown that the pressures on a cylinder with ; 


vortices symmetrically situated downstream (as is shown in Oa ; n ~ aor . : 
Fig. 1) depend upon the position and the motion of the vortices * at r Ot waist a 
relative to the center of the cylinder. The complex potential a 
function for this flow is 2elv + 2eP be (u sin 2y — vcos 2y) (6) 


W (¢) U (« T 2 c 2 where u and v are the velocities of the vortices in the direction of 
— éand 7, respectively. The ¥ force is zero from symmetry. 


The first term in Eq. (6) is identical to the result obtained by 
Munk? and the last two terms are Michelson’s results. This 
more general expression indicates, however, that the force is 
aie altered if the cylinder is growing in the presence of vortices and /or 

+ in the vortex strength is changing with time 
"Y = x + Ly, instantaneous position of the vortex with 
strength F. REFERENCES 


The Og 2 evi ar oo > eYE ‘ > : yrati > : : 
rhe force on the cylinder can be evaluated by integrating the ; nikitidieies 6 aie wie den of Valtent Diliel « Che Cee 
pressure around the boundary from Naval Ordnance Test Station, Inyokern, Calif., TM 377, 1951 

2 Munk, Max M., The Aerodynamic Forces on Airship Hulls, NACA Rep 


X¥ -iY= f — p(dn + idé) : 184, 1924. 
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